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ABSTRACT 

Consider  a general  Two-Point  Boundary  Value  Problem  (TPBVP) : 
y'  (t)  = f (t,y)  ^ t 1 

B^y(a)  + B2y{b)  = w 

where  frR"^  ^ ^ R*^,  f e C^,  B^  and  B^  are  nxn  matrices  and  we  R*^. 

It  is  shown  how  one  can  bound  a posteriori  the  error  made  in  the  numerical 
solution  of  the  TPBVP.  The  error  bounds  obtained  are  rigorous  and  include  the 
truncation  and  the  roundoff  error.  In  addition,  the  confutations  establish 
the  existence  of  solutions  to  the  TPBVP.  ' 

Numerical  schemes  are  developed  for  the  case  where  f(t,y)  is  a polynomial 
in  t and  y . Examples  are  given  of  conf utational  existence  proofs  for 
problems  where  analytical  existence  proofs  are  not  known. 
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SIGNIFICANCE  AND  EXPLANATION 


Consider  a general  two-point  boundary  value  problem 


y'  (t)  = f(t,y (t) 


B^y  (a)  + B^y (b)  = w 


a < t < b 


where  y(t)  and  f(t,y(t))  are  n dimensional  vector-valued  functions, 

Bj^  and  B^  are  nxn  matrices  and  w is  n vector. 

Many  problems  in  physics,  chemistiry  and  engineering  can  be  put  into  the 
above  form.  In  most  cases,  it  is  impossible  to  write  down  explicit  formulas 
for  solutions  of  two-point  boundary  value  problems.  Moreover,  many  times  it 
is  impossible  to  assert  the  existence  of  a solution  or  to  know  whether  the 
equations  have  one,  many  or  an  infinite  number  of  solutions. 

Although  there  are  many  numerical  methods  that  enable  one  to  compute 
numerical  approximations  to  solutions  of  equation  (1),  it  is,  in  general 
impossible  to  tell  a priori  how  close  the  numerical  solution  is  to  a true 
solution . 


The  literature  holds  examples  of  two-point  boundairy  value  problems  that 
do  not  have  a solution  but  the  equations  arising  from  the  numerical  approxima- 
tion scheme  do  have  solutions.  Also,  there  are  examples  of  numerical  solutions 
that  at  first  were  believed  to  be  good  approximations  but  later  were  proved 
to  be  wrong. 

We  developed  a numerical  scheme  that  enables  one  to  take  a numerical 
solution  to  equation  (1)  and  with  the  aid  of  that  solution  to  compute,  a 
posteriori,  guaranteed  error  bounds.  That  is,  one  can  compute  how  far  the 
numerical  solution  is  from  the  true  solution. 

The  method  described  in  this  report  is  applicable  in  the  cases  where 
f(t,y)  is  a polynomial  in  t and  y . Although  this  is  a serious  restriction, 
many  problems  that  arise  in  applications  do  fall  into  that  category.  We  also 
give  some  indication  of  how  one  would  extend  our  procedure  to  more  general 
types  of  two-point  boundary  value  problems,  that  is,  cases  where  f(t,y) 
is  not  a polynomial  in  t and  y . 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 
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CHAPTER  1 


A POSTERIORI  ERROR  BOUNDS 
FOR  TWO-POINT  BOUNDARY  VALUE  PROBLEMS 


1.1 


Introduction 

Consider  a general  Two-Point  Boundary  Value  Problem  (TPBVP) 
y' (t)  = f (t,y(t) ) 0 1 t 1 1 


(1) 


B^y(O)  + B^yd)  = w 

, ..n+1  n - 2 , „ „ n 

where  f:R  ^ R , f e C and  B^,  B^  are  nxn  matrices,  w e R . 

Unlike  Initial  Value  Problems  where  the  theory  assures  us  that 

for  a large  class  of  equations  there  is  a unique  local  solution, 

TPBVP  can  have  one  or  many  or  no  solutions  at  all.  Unless  one  makes 
very  strong  assunptions  on  the  function  f in  equation  (1)  the  ques- 
tion whether  a solution  does  exist  or  whether  it  is  unique  is  hard  to 
answer.  Many  times  numerical  solutions  to  TPBVP  are  computed  without 
establishing  the  existence  of  a solution.  The  error  analysis  of 
numerical  methods  for  solving  TPBVP  assures  us  that  if  a solution 
does  exist  cind  is  locally  unique  and  the  discretization  parameters  are 
small  enough,  there  is  a locally  unique  numerical  solution  which  is 
close  to  the  true  solution  (see  Keller  [22],  de  Boor  and  Swartz  [2], 

Krasnoselskii  et  al  [24]).  The  above  theory  is  unsatisfactory  for 

two  reasons;  a)  It  is  hard  if  not  in^jossible  to  prove  by  analytic 
methods  that  there  is  a solution,  b)  In  general  it  is  impossible  to 
know  how  small  to  make  the  discretization  parameters. 

Several  numerical  methods  for  estimating  the  error  made  in  the 
numerical  solution  of  differential  equations  are  known.  The  most 
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famous  one  is  Richardson  Extrapolation  to  the  Limit  approach:  If 

the  equations  and  the  solution  are  smooth  enough,  one  can  show  that 

the  error  admits  an  expansion  of  the  form 

k 

y - yU  ) = I h^H.  (X  ) + OCh'^  . 
n n in 

D=P 

Based  on  this  expansion  one  can  get  asymptotic  error  estimates 
and  one  can  iitprove  the  solution.  In  [47],  Pereyra  gives  a survey  of 
such  methods.  Honrici(see  [161)  points  out  that  one  can  solve  numer- 
ically the  differential  equation  satisfied  by  the  dominant  error  term 
in  order  to  get  a good  estimate  of  the  error.  Some  methods  like  the 
Runga-Kutta-Fehlberg  methods  have  error  estimates  built  in.  In  [46] , 
Zadunaiski  describes  another  type  of  method.  He  interpolates  the 
numerical  solution  by  a local  polynomial.  Then,  using  that  poly- 
nomial, he  constructs  a pseudo-problem  for  which  the  solution  is 
known.  He  estimates  the  error  in  the  numerical  solution  by  computing 
the  error  made  by  the  integration  process  in  solving  the  pseudo- 
problem. More  recently  Babuska  and  Rheinboldt  [ 1 ] analyzed  some 
aposteriori  error  estimates  for  TPBVP  of  elliptic  type  solved  by  the 
finite  element  method. 

Numerical  estimates  usually  work  very  well.  However,  these  are 
only  estimates.  Numerical  estimates  can  always  fail.  It  was  pointed 
out  by  Lyness  and  Kaganove  in  their  paper  "Comments  on  the  nature  of 
Automatic  Quadrature  Routines"  [25],  that  numerical  schemes  which 
estimate  the  error  by  considering  function  values  at  finitely  many 
points  inevitably  may  fail,  even  in  the  simplest  case  y' (t)  = f(t), 
that  is  integration. 


i 


Therefore,  it  is  desirable  to  construct  numerical  schemes  that 


will  enable  one  to  take  a numerical  solution  and  bound  aposteriori 
the  error  in  that  solution.  By  doing  so,  one  also  proves  existence 
of  a solution.  We  would  like  to  stress  that  we  are  not  concerned 
with  error  estimates.  We  are  interested  in  computing  guaranteed  error 
bounds.  Error  bounds  are  more  expensive  than  error  estimates,  but  in 
many  cases  the  ability  to  talk  about  the  accuracy  of  the  numerical 
solution  with  certainty  is  inportant. 

Many  times  the  existence  of  a solution  to  a single  second  order 
Boundary-Value  problem  with  separated  boundary  conditions  can  be 
established  by  the  "Shooting  Method":  One  can  find  two  values  of  the 
missing  initial  condition  such  that  the  solution  of  the  initial 
value  problem  with  one  initial  value  will  "hit"  above  the  end 
condition  and  the  other  below  it.  Since  the  solution  of  the  differ- 
ential equation  is  continuous  as  a function  of  the  initial  conditions 
(see  Coddington  and  Levinson  [ 8 ]),  there  must  be  a value  of  the  mis- 
sing initial  conditions  that  will  make  the  solution  "hit"  the  right- 
hand  side  exactly.  However,  the  above  idea  cannot  be  used  to  analyze 
more  general  problems. 

A very  general  approach  to  the  question  of  existence  and  apos- 
teriori error  analysis  was  suggested  and  analyzed  by  L. V.  Kantorovich. 
He  generalized  Newton's  method  for  solving  nonlinear  equations  to  a 
general  method  of  solving  nonlinear  operator  equations  in  Banach 
spaces. 
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In  order  to  describe  Newton's  method  in  Banach  space  and 
Kantorovich's  theorem  we  first  describe  the  following  concepts'. 

Let  X and  Y be  Banach  spaces  and  L be  a bounded  linear  operator 

from  X to  Y.  The  induced  operator  norm  of  L is  defined  by 

II  l||  = sup  II  Lx||y  . 

II  X 11^=1 

Let  B(X,Y)  be  the  set  of  all  bounded  linear  operators  from  X to 

Y.  If  we  define  addition  and  scalar  multiplication  in  B(X,Y)  by 

V=(W+Z)‘*  Vx=Wx  + Zx  VxcX. 

V = aZ  «■  Vx  = aZx  V x e X 

then  it  is  easy  to  see  that  V is  a bounded  linear  operator  from  X 

to  Y and  that  B(X,Y)  is  a vector  space.  It  is  well  known  that 

B(X,Y)  is  itself  a Banach  space.  (see  [20]). 

Let  P be  an  operator  from  an  open  set  D 5^  X into  Y . If 

there  is  a bounded  linear  operator  L e B(X,Y)  such  that 

II  P(x  +u)  - P(x  ) - Lu|L 
lim  = 0, 

Wlx 

then  P is  said  to  be  differentiable  at  x^  and  the  bounded  linear 
operator  P'  (x^)  = L is  called  the  Frechet  derivative  of  P at  x^. 
If  P is  differentiable  at  every  point  of  the  open  set  D then  the 
mapping  P'  : x^  P'  (x^)  is  an  operator  from  D c X to  B(X,Y) 
and  therefore  one  can  talk  about  the  derivative  of  such  an  operator. 
That  derivative  is  called  the  second  Frechet  derivative  of  P,  at  x^, 
and  is  denoted  by  P" (x^) . 

^For  a full  treatment  see  Kantorovich  and  Akilov  [20]  cind  also  Rail 
[32]. 
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Note  that  P"(Xq)  e B(X,B(X,Y))  and  for  u,v  £ X, 
P"  (x^)  [u]  e B(x,y)  and  (P"  (x^)  [u] ) (v]  e Y . 

The  induced  operator  norm  of  P"{Xq)  is: 


P"(Xg)ll  = su: 


“1= 


sup  ||{P"(X  ) [u])  [V]  |l 

1 llv||j^=l  "" 


Example : 


Equation  (1)  can  be  put  in  the  above  framework.  Let  X^  be  the 
space  (C^[0,1])*^  and  X^  = (C  [0,1])’’  x that  is: 

, I 

X 


e C^[0,1],  i = l,...,n} 
X,  . ([^"]  « C(0,11 


, P £ R } . 


A norm  on  is  (for  example) 

II  f II  = max  [max(  sup  |f.(t)|,  sup  |f ! (t)  | )]  . 
i lliln  0£t^l  0£t<l  ’■ 

A norm  on  X^  is  (for  example) 

max(  max  [ sup  (|f.(t)|)],  max  Ip.j)  . 
l^i^n  O^t^l  ^ ^ 

It  is  not  hard  to  show  that  with  the  above  definitions  X^^  and  X^ 


are  Banach  spaces. 


Let  P : Xj^  -*■  X2  be  defined  by: 


P(x) (t)  = 


x'  (t)  - f (t,x(t)) 


B^x(O)  + B^xd)  - w 


0 < t < 1 . 


Clearly  solving  for  P(x)  =0  is  equivalent  to  solving  equation  (l) 
P' (Xq)  is  the  linear  operator  defined  by 


P' (XQ)v(t)= 


v' (t)  - f (t,x  (t) ) v(t) 
X 0 

B^v(0)+  B^vd) 


where  f (t,x„(t))  is  the  matrix 
X 0 


3 f 


3x.  y.  • 1 ' 

3 ' id=l 


0 < t < 1 
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< 


The  second  Frechet  derivative  is  given  by: 


P" (Xq)  [v,u]  (t) 


(f  (t,x^(t) )v(t) )u(t) 

XX  0 


0 < t < 1 


where  f (t,x„(t')  is  the  tensor 
XX  0 


3 f.(t,x^(t)) 
1 0 

3 X , 3 X 
D k 


Note 


i, j,k=l 

y-x|l  is 


that  the  norm  in  involves  the  first  derivative  and 

small  only  if  sup  |y(t)-x(t)|  and  sup  | y' (t) -x' (t) | are  small. 
0^t<l 

We  are  now  ready  to  describe  Newton's  method  in  Banach  space: 

Recall  the  regular  Newton's  method:  One  is  trying  to  find  x^  such 

that  f(x  ) = 0.  Assume  that  we  have  x„  which  is  an  approximation 

0 

to  Xj^ . We  try  to  find  x^  which  is  a better  approximation  to  x^  . 
We  write, 

0 = f (x  ) = f(x  +Ax)  = fix  ) + f (x  )Ax  + Q{Ax) 
u u u 

If  one  neglect  the  term  Q(Ax)  and  solve  for  Ax  one  gets 
Ax=-f  (Xq) /f ' (Xq).  Now  if  x^  = x^  + Ax  then  x^  = x^-f  (x^) /f ' (x^^)  , or 
in  general  x^^j^=  x^^ -f  (Xj^)/f ' (x^^)  which  is  Newton ' s iteration  scheme. 

Following  the  same  reasoning  one  can  generalize  the  above  proce- 
dure to  operator  equations  in  Banach  space.  Let  x^  be  an  approxi- 
mation to  x^ . Assume  that  P is  twice  continuously  differentiable. 
We  write 


0 = P(x.)  = P(x  +Ax)  = P(x^)  + P'(x„)Ax  + Q(Ax)  . 
"0  0 u 

By  do^'inition  of  Frechet  derivative,  P'  (x  ),  lim  iL 2 (A?  1.11 

° II  AxlUo  II  Ax  II 


= 0 


So  if  II  Ax  II  is  sufficiently  small  one  can  neglect  the  operator 
Q(Ax)  and  solve  for  Ax, 

Ax  = - [P'  (x^)  ] (Xq)  , 

assuming  of  course  that  [P' (x^]  ^ exists.  Again  the  general  scheme 
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is  X = X.  - [P'(x  )]  P(x,). 

k+1  k k k 

Analyzing  the  above  scheme  Kantorovich  proved  the  following  re- 
markable theorem: 

Theorem:  (Newton's  method  in  Banach  space). 

Let  X and  Y be  Banach  spaces  and  let  P:D  ^ X Y (D  open 

2 

set)  be  a nonlinear  operator.  Assume  that  P e C (X,Y) . Let 
Xq  t X and  assume  the  following: 

a)  [P' (Xq) ] ^ exists 

b)  II  Xj^-XqII  ^ n where  “ [P'x^))  ^P(Xq) 

II  [P’  (Xq)  ]"^|1  < B 

d)  ||p"(x)||  < K for  ||x-Xq|(  <r 

e)  h = n-  B • K < |- 

1-  A-2h 

f)  ^1^0  = h ^ 

then  P(x)  = 0 has  a unique  solution  x^  in  the  ball  ||  x-x^H  ^ r^ 

and  the  iteration  scheme  x ,=  x - [P' (x, ) ] ^P(x, ) will  converge 

K+ 1 K -K  K 


The  above  theorem  enables  one  to  prove  existence  by  first  com- 


puting an  approximate  solution,  then  by  computing  the  relevant  con- 
stants one  can  establish  the  existence  of  a solution  nearby.  In 
addition,  estimate  (f)  gives  bounds  on  the  difference  between  the 
computed  solution  and  the  exact  solution.  Although  the  theory  has 
been  known  for  more  than  twenty  years  it  has  been  used  very  little 
in  order  to  establish  the  existence  of  solutions  to  continuous  prob- 
lems like  TPBVP.  After  a moment  of  thought  one  realizes  that  the 
actual  application  of  the  above  abstract  results  to  TPBVP  is  not 
trivial.  A careful  formulation  of  the  problem  has  to  be  made. 
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In  order  to  compute  a rigorou  bound  on  the  error  one  has  to  be 
able  to: 

a)  Compute  a reasonable  bound  on  |1  [F'  1 ^11 

b)  Compute  a good  bound  on  |1  x^-x^H  . 

c)  Compute  a bound  on  sup  ||  F"(x  ) ||  . 

II  x-x„||  < . “ 

First,  one  has  to  choose  the  spaces  X and  Y . A natural  choice  is 

and  X^  of  the  example.  However,  since  the  norm  in  X^  depends 
on  the  first  derivative  and  one  needs  to  make  h = |1  small, 

the  above  choice  implies  that  to  be  a good  approximation 

to  x^(t).  This  is  a very  stringent  requirement.  Therefore,  it  is 
desirable  to  find  a formulation  and  a space  X where  the  norm  of 
X depends  on  function  values  only.  Even  in  that  case  bounding  the 
norm  requires  bounding  the  range  of  functions;  hardly  a trivial  mat- 
ter. (See  Moore  [29],  and  the  references  there.)  Second,  computing 
a bound  on  j]  F' (x^)  ^|1  implies  bounding  of  norm  of  the  Green's 
function  of  a first  order,  linear,  two-point  boundary  value  problem. 
Theoretical  bounds  for  such  a problem  are  usually  impossible  to  ob- 
tain so  one  has  to  devise  a conputational  scheme  for  getting  such  a 
bound.  Third,  one  has  to  confute  n and  K . This  computation  takes 
bounding  the  range  of  one  and  two  dimensional  functions  on  a finite 
domain.  This  is  not  an  easy  task.  (See  Moore  [29]  and  the  refer- 
ences there.)  All  of  the  above  gets  even  more  complicated  by  the 
fact  that  computers  compute  with  only  a finite  number  of  digits  and 
every  arithmetic  operation  introduces  some  round  off  errors.  If  one 
wishes  to  prove  existence  of  a solution  one  has  to  take  that  point 
into  account. 
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The  residuals  will  be  very  small  functions  that  are  the  differ- 
ence between  almost  equal  but  not  small  functions.  Therefore^  it  is 
impossible  to  decide  apriori  how  many  digits  one  needs  to  use  in 
order  to  be  able  to  neglect  the  effect  of  round-off  error. 

Several  authors  have  suggested  various  formulations  and  compu- 
tational schemes  that  enable  one  to  compute  the  constants  of 
Kantorovich's  theorem.  However,  so  far,  the  methods  that  are  sug- 
gested in  the  literature  have  some  difficulties.  Talbot  [37]  has 
constructed  a method  that  enables  one  to  use  a modified  version  of 
Kantorovich's  theorem  to  bound  the  error  for  single  second  order 
boundary  value  problem  y"  = f(t,y)  where  the  right  hand  side  does 
not  depend  on  y' . In  order  to  take  the  truncation  and  the  round- 
off error  into  account  Talbot  uses  Interval  Analysis  techniques 
(see  Moore  [30]).  He  computes  an  interval  valued  function  that 
contains  the  exact  Green's  function  in  its  range.  By  computing  a 
bound  on  the  range  of  the  interval  function,  he  is  able  to  compute 
a bound  on  ||  F'Cx^)  ^||  . He  also  uses  interval  techniques  to  com- 
pute bounds  on  ||  ||  F"||  . His  approach  enables  him  to 

compute  guaranteed  global  bounds  on  the  error.  The  major  disad- 
vantage of  his  method  is  that  it  requires  large  memory  storage  space. 
Also,  as  pointed  out  by  Talbot  himself,  it  is  desirable  to  solve  a 
more  general  class  of  problems.  McCarthy  and  Tapia  [26]  have  sug- 
gested a different  approach.  They  convert  a single,  polynomial, 
second  order,  two-point  boundary  value  problem  into  an  equivalent 
Volterra  equation  and  use  a modified  version  of  Newton's  method,  due 
to  Tapia,  to  confute  error  bounds  on  the  approximate  solution. 
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Their  method  also  suffers  from  some  difficulties.  The  bound  on 

II  F'  (Xq)  ^11  is  given  jn  terms  of  6^  where  L is  a bound  on  the 

norm  of  a certain  function.  This  bound  is  very  pessimistic.  For 
L 8 

L = 20 , e > 10  . Such  a large  bound  makes  the  possibility  of  get- 
ting realistic  arror  bounds,  or  error  bounds  at  alJ,  rather  small. 
Also  the  autho'-s  do  not  consider  the  effects  of  round-off  error  at 
all. 

Roberts  and  Shipman  [34]  have  suggested  to  convert  the  problem 
of  two-point  boundary  value  problem  to  a set  of  algebraic  equations 
by  looking  for  the  missing  initial  conditions.  Then,  they  solve 
that  set  of  equations  by  the  "shooting  method".  However,  they  do 
not  consider  the  error  made  by  the  process  of  solving  the  initial 
value  problem  and  therefore,  their  approach  cannot  be  used  to  get 
rigorous  error  boiands. 

In  [12]  Cruickshank  and  Wright  suggested  a method  for  bounding 
11  F'(x  ) ^11  for  2mth  order  single  TPBVP.  They  use  the  Green's 

<J 

function  of  the  2mth  derivative  plus  boundary  conditions  to  convert 
the  TPBVP  into  an  equivalent  Fredholm  integral  equation  of  the 
second  kind.  They  use  Kantorovich's  Theory  of  Approximation  Methods 
(see  [20],  Ch.  XIV),  to  relate  the  norm  of  (I-k)  ^ to  the  inverse 
of  the  matrix  that  one  gets  by  using  polynomial  collocation  method 
to  solve  the  TPBVP.  Their  approach  can  also  be  used  for  systems  of 
equations  and  piecewise  polynomial  collocation  methods  as  well. 
However,  their  approach  can  lead  to  the  inversion  of  very  large 
matrices.  Although  the  matrices  themselves  are  band  or  block  band 
matrices,  their  inverse  is  full.  Also  it  is  not  always  possible  to 
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find  an  explicit  form  of  the  Green's  function  that  will  enable  one 


to  convert  the  problem  to  an  equivalent  Fredholm  integral  equation. 

Our  method  is  to  convert  equation  (1)  into  an  equivalent  Vol terra 
equation.  Then  we  derive  an  explicit  formula  for  the  Green’s  func- 
tion of  that  Volterra  equation.  Using  the  formula  we  construct  a 
numerical  scheme  to  compute  bounds  on  ||  F' (x^)  ""  )|  and  || 
for  function  f(t,y)  which  are  polynomials  in  t and  y . We  are 
using  only  spaces  of  continuous  functions. 

This  has  two  advantages; 

a)  We  can  use  local  polynomial  interpolation  in  order  to  produce 

continuous  functions  as  opposed  to  global  piecewise  polynomial 

functions . 

b)  Only  the  values  of  the  computed  solution  have  to  be  a good 
approximation  to  the  true  solution  and  not  to  the  derivatives. 
Getting  a good  approximation  to  the  solution  and  its  derivative  is 
of  course  harder  than  obtaining  a good  approximation  only  to  the 
solution  itself.  The  requirement  that  f(t,y)  in  equation  (l)be  a 

polynomial  in  t and  y is  very  restrictive.  However,  a large 

I 

number  of  problems  that  arise  in  applications  (most  notably  one 
dimensional  cases  of  Naviar  Stokes  equations) , fall  into  that  cate- 
gory. 

We  also  would  like  to  note  that  our  method  can  be  made  to  work 
for  functions  f(t,y)  which  are  factorable  functions  of  t and  y, 
f(t,y)  can  be  approximated  by  Taylor  series  approximation  and  tlie 
■ error  made  by  that  approximation  can  be  accounted  for.  However,  we 

do  not  treat  that  case  and  we  leave  it  for  further  study. 
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I 


We  conclude  the  introduction  with  an  example  of  TPBVP  that  can 
be  treated  by  our  method.  The  example  looks  simple  but  is  far 
from  being  trivial. 

Consider  the  system  of  ordinary  differential  equations  des- 
cribing the  steady  flow  of  a fluid  contained  between  two  parallel, 
infinite  plane  disks  which  are  rotating  about  a common  axis 
h^'^  + hh'"  + gg'  =0 
g"  + hg'  - h’g  = 0 
h(0)  = h' (0)  = h(l)  = h’ (1)  =0 
g(0)  = , g(l)  = . 

The  above  problems  have  attracted  considerable  attention. 

In  [27],  McLeod  gives  a list  of  more  than  twenty  references  about 
work  concerning  the  above  problem.  Although  the  physically  inter- 
esting problem  is  when  | | and  | | are  very  large,  existence 
proofs  are  known  for  a very  limited  set  of  values  of  (see 

Elcrat  [14]).  Moreover  extensive  computations  have  been  carried 
out  and  multiple  solutions  had  been  observed  (see  Wilson  and  Schryer 
[4ll  and  the  references  there).  However,  the  known  theory  cannot 
assert  the  existence  of  such  multiple  solutions  nor  could  it  rule 
them  out. 

One  "nice'  feature  of  the  above  equation  is  that  it  is  q ladratic 
in  the  unknown  functions.  Therefore,  F'  is  a constant  tens.-r. 

Using  our  method  we  prove  the  existence  of  a solution  for  boundary 
conditions  outside  the  range  for  which  analytical  proof  is  known. 
Actually,  Kantorovich's  theorem  guarantees  the  existence  of  a solu- 
tion for  boundary  values  in  some  ball  around  the  original  boundary 
values. 
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1 . 2 The  Theory : 


In  order  to  apply  Kantorovich's  Theorem  to  the  problem 
y'  = f(t,y) 


{ 


0 < t < 1 


B^y(O)  + B^yd)  = w 


(1) 


we  use  the  following  formulation: 

Let  C[0,1]*^  denote  the  space  of  continuous  functions  from  [0,1]  to 

R . Let  r e We  use  |r|  to  denote  max  |r.  |.  For  f £ 

1 ^i  ^ n ^ 

we  define  the  norm  of  f by; 

II  f II  = sup  I f(t)|  = sup  max  (|f.(t)|). 


0 < t < 1 


0<t<l  l<i<n 


We  denote  by  Cq[0,1]  the  space 

Cq[0,1]"  = {g  € cro,!]*^  I g(0)  = 0}  . 
Clearly  C[0,1]^  and  Cq[0,1]'^  with  the  norm  1(*|1 


are  Banach  spaces, 
n 


Let  A be  n^n  matrix.  We  use  |a|  to  denote  max  ( ^ ]a. ,|). 

l^i<n  j=l 

Let  W(t)  be  a nxn  matrix  valued  function  on  [0,1].  We  define  the 
norm  of  W by 


I W II  = sup  I W(t)  I . 
0 < t<  1 


Let  f be  a twice  continuously  differentiable  function  from  R*^^^  to 
r"^.  Define  F:C[0,1]'"  ->  €^[0,1]"  r''  by 


f u(t)  - u(0)  - /^f (s,u(s) )ds 
Fu(t)  = ( 0 

l.B^u(0)  + B^ud)  - w 


(2) 


Clearly  solving  for  Fu  = 0 is  equivalent  to  solving  equation  (1)  . 
f”  (Xg)  : the  Frechet  derivative  of  F at  x^  is  given  by: 

v(t)  - v(0)  - f (s,x  (s))v(s)ds 

_ X u 


F’  (Xg)  [V]  (t)  = 


(3) 


B^v(O)  + B^vd) 
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n 


where  f (srx)  is  the  matrix 

X 


3x.  I . . , 
: J i.:=i 


rf  n c 


n = ^ then  V = F*  (x  ) [n]  is  the  solution  of 

V p 0 


v(t)  - v(0)  - f (s,x  (s))v(s)ds  = r(t) 
0 ^ ° 


B^v{0)  + V (1)  = p 


F"  (Xq)  the  second  Frechet  derivative  is  given  by: 


F"  (X|^)  [v 


,w] (t)  = < 0 

L ° 


(f  (s,x  (s) )u(s) )w(s)  ds 

XX  0 


where  f (t,x)  is  the  tensor 

XX 


j,k  — l,...,n 


Let  Xq  e C[0,1]  be  fixed  and  denote  f^{s,x^{s))  by  A(s).  A(s) 
is  a nxn  matrix  valued  function. 

Let  Y{t)  be  the  n^n  matrix  that  solves  the  following  initial 


value  problem: 


Y’ (t)  = A(t) Y(t) 


Y(0)  = I 


0 < t < 1 


We  are  now  ready  to  derive  a formula  for  F' (x^) 
Theorem: 

Let  R = B^  + B^  Y(l) 

F' (Xg)~^  exists  if  and  only  if  R is  nonsingular  and 


" y(t)R'\/^ 


Y (s) A{s) r (s)ds 


+ Y (t)  (r"^Bj^-I) Y (s)A(s)r(s) 


- Y(t)R  (B^  r(l)  + p)  + r(t) 


Proof : 


Since  A(t)  is  continuous  on  [0,1]  it  is  well  known  (see 


Coddington  and  Levinson  [8])  that  Y(s)  is  nonsingular  for  all 


I 


s ■ [0,1).  So  equation  (7)  is  well  defined  provided  R ^ exists. 

Let  w(t)  = Y(t)  Y“^(s)A(s)r(s)ds  + r(t),  then,  w(0)  = 0 and 
0 ^ 

w(t)  - w(0)  - f A(s)w(s)ds  - r(t)  = 

Y(t)  y“^ (s)A(s)r(s)ds  - A(s)Y(s)  y“^ (u) A (u) r (u) du  ds 

0 0 0 

- A(s)r (s)ds  = 

0 

(integrating  the  middle  terms  by  parts) 

= Y(t)/^  Y"^(s)A(s)r(s)ds  - Y(t)  Y-^(u)A(u)r(u)du 
0 0 

+ y(s)Y  ^ (s)A(s)r(s)ds  - A(s)r(s)ds  = 0 
0 0 

therefore 


w(t) 


w(0)  - A(s)w(s)ds  = r(t)  . 


0 


Since  Y(t)  - y(0)  - A(s)Y(s)ds  = 0,  for  any  a e 

0 

v(t)  = w(t)  + Y(t)a  is  a solution  to: 

v(t)  - v(0)  - A(s)v(s)ds  = r(t).  (8) 

0 

Clearly  any  solution  to  (8)  is  of  the  form  w(t)  + Y(t)a  for  some 

u £ r"^.  If  V is  a solution  of  equation  (4)  then 
rt  -1 


V 


(t)  = Y(t)  / Y * {s)A(s)r(s)ds  + Y(t)a  + r(t) 
0 


n 


for  some  a £ R 


B^v(O)  + B^vd)  = p 


The  condition 
implies  that 


B a + B [y(l){  Y ^ (s)A  (s)  r (s)ds  + a}+r(l)]=  P 
^ 0 

hence 

(B  + B Y(l))a  = p-B  [Y(l)  y"^(s)A(s)r (s)ds  + r(l)l 

0 

hence 

Rf«  = p-B  [Y(l)  Y~^(s)A(s)r(s)ds  + r(l)] 

0 


(9) 


Therefore  if  R is  nonsingular 

a = R ^ {p  - B [Y(l)  Y ^ (s) A (s) r (s) ds  + r(l)]) 

0 


v(t)  = Y(t){/  Y (s)A  (s)  r (s)  ds  + R [p-B  Y(l)/  Y ( s)  A ( s ) r ( s)  ds 
0 0 

- B^r  (1)  ] } + r(t) 


= Y(t)R  ^ {(B^+B^Y (1) ) /^Y  ^ (s)A(s) r (s)ds 


B Y(l)/  Y (s)A(s)r (s)ds  + p-B  r(l)}  + r(t) 


= Y(t)R~^B  /S"^(s)A(s)r(s)ds  - Y(t)R"‘^B  Y(l)/S"'^(s)A(s)r(s)ds 
^0  t 


+ Y(t)R  -"(p-B^rd))  + r(t) 


= Y(t)R  ^B^/^Y  ^ (s) A (s) r (s) ds 
■^0 

-Y(t)R  ^(R-B^)/^Y  ^ (s)  A (s)  r (s)  ds  + Y(t)R  ^(p-B^rd))  + r(t) 


= Y(t)R  ^B^/^Y  ^ (s)A(s)  r(s)ds  + Y(t)  (R  ■'■B^-I)/'’'Y  (s)  A (s)  r (s)  ds 


+ Y(t)R  "(p-B^rd))  + r(t)  . 

Conversely  by  virtue  of  the  above  computation,  equation  (9)  has  a 


unique  solution  only  if  R is  nonsingular. 


Q.E.D. 


Remarks : 


1)  If  we  replace  Y ^(s)A(s)  by  -(Y  ^(s))'  , we  can  rewrite  equa- 
tion (7)  as  : 

F'  = - Y(t)R~^Bj^(Y"^(s))  •r(s)ds  - 

P • ^0 

- Y(t)  (r"^B  -I)/^(y“^(s)  ) 'r(s)ds  - Y(t)R"^(B  rd)+p)  + r(t).dl) 
^ t ^ 

2)  Using  the  above  equation  we  get 

II  F'(x  )”^|1  < ||y||  {(  /^|  (Y"^(s))'|ds)max(|R‘^B  1,|r"^B  -i|)  + |r"^B  I 
° ■ 0 ^ 

+ |R~^i } + 1.  (12) 


i 


1.3  Bounding  the  Constants  of  Kantorovich's  Theorem: 


In  order  to  construct  an  effective  numerical  scheme  we  assume 
that  the  function  f(t,y)  in  equation  (1)  is  a polynomial  in  t 
and  y . 

We  would  like  to  note  that  all  functions  involved  can  be  ap- 
proximated by  local  Taylor  series  expansions  and  therefore  the  pro- 
cedures we  are  about  to  describe  can  be  used  to  compute  the  desired  1 

bounds  even  if  f is  not  a polynomial.  The  difference  between  the 
exact  equation  and  the  local  piecewise  polynomial  approximation  can 
be  accounted  for.  In  the  present  work  we  do  not  investigate  this 
possibility,  and  we  leave  it  for  further  study.  In  any  case  as  was 
pointed  out  before,  many  problems  that  arise  in  applications  are  poly- 
nomials in  t and  y . 

In  order  to  take  into  account  the  round-off  error  we  are  using 
interval  arithmetic.  The  use  of  interval  arithmetic  at  the  present 
time  is  very  expensive  since  one  has  to  use  simulated  floating  point 
arithmetic.  The  ratio  of  speed  between  the  hardware  floating  point 
arithmetic  and  the  interval  arithmetic  we  are  using  is  about  1:1000. 

Although  the  scheme  we  are  proposing  does  not  require  a large 
volume  of  computation,  the  use  of  interval  arithmetic  makes  the  method 
expensive  to  use.  If  one  is  not  interested  in  actually  proving  the 
existence  of  a solution  and  one  does  not  wish  to  take  round-off  error 
into  account  the  use  of  double  precision  will  suffice.  We  would  like 
to  point  out  that  in  computing  ||  F"(Xq)  ^||  and  ||  we  use 

interval  arithmetic  only  in  order  to  take  the  round-off  error  into 
account.  Therefore,  we  do  not  get  the  pessimistic  bounds  that  the 
use  of  interval  arithmetic  can  lead  to. 
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1.3.1  Bounding  ||  F"  (x)  ||. 


If  the  function  f(t,y)  in  equation  (1)  is  a polynomial  in  t 

and  y and  we  choose  our  initial  guess  ^ piecewise 

polynomial  function,  then  f (t,x„(t))  will  be  a tensor  of  order  3 

XX  0 

whose  components  are  piecewise  polynomial  functions. 

One  way  tc  compute  K (a  bound  on  sup  ( l|f  (t,x(t))  jj) 

II  X-Xgll  ir 

is  to  bound  the  norm  of  F (t,x  (t)-I)  where  F is  the  interval 

XX  0 XX 

extension  of  f^^  (see  Moore  [30])  and  I is  the  interval  [-r,r]. 

[F  (t,x  (t)-I)].  , is  a piecewise  polynomial  intejrval  function. 

XX  U 1 / ^ / X 

In  section  (1.4)  we  describe  how  one  can  bound  the  norm  of  such 
functions.  Another  way  to  compute  K is  to  bound  (|  f^^  (t,X|^(t))|| 
and  then  to  coirpute  a crude  bound  on  sup  |1  f (t,x(t))|l  . 

Since  each  component  of  f (t,x)  is  a polynomial  in  t and  x 

XXX 

it  is  not  hard  to  compute  a crude  bound  on  H f ^^^(t,x(t))  j]  . 


1.3.2  Bounding  ||  F' (x^)  ^||  . 

Let  A(t)  denote  f^(t,XQ(t)).  Let  Y(t)  be  the  matrix  solu- 
tion of  the  initial  value  problem 


Y' (t)  = A(t)  • Y(t) 


0 < t < 1 


Y(0)  = I 

It  is  well  Itnown  that  Y ^ (t)  exists  and  is  the  solution  of: 
Z*  (t)  = -Z(t)A(t) 

Z(0)  = I . 

Let  V(t)  and  W(t)  be  matrix  valued  functions  approximating 
and  Y ^(t)  respectively.  Assume  V(0)  = W(0)  =1  . 

Let  E^(t)  = Y(t)  - V(t),  E2(t)  = Y ^(t)  - W(t).  satisfies 


Y(t) 


the 
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equation 

E^(t)  = A(s)E^(s)ds  + X (t) 

0 . 

where  X (t)  = I - V(t)  + / A(s)V(s)ds  also 

0 

E (t)  = E (s)A(s)ds  + X (t) 

0 

where  X (t)  = I-W(t)-/^  W(s)A(s)ds.  Therefore 
^ 0 

!E^(t)|  < |x^(t)l  + |a(s) I 1e^(s) Ids 

In  order  to  bound  the  norm  of  E.  in  terms  of 

1 


(13a) 


(13b) 


i = 1,2. 

{ , II  and  II  A | 


one  can  use  the  following  well  known  lemma  (See  [8]  Chap.  1). 
Lemma  (Gronwall  inequality) . 

Let  'firVrX  be  real  piecewise  continuous  functions  on  a real 

interval  I = [a,b]  , 0 on  I and  suppose  that  Vt  € I 

■,p(t)  ^ X (s)'fi  (s)ds . 

a 

Then 

>^(t)  <_  i|/(t)  + x(s)t|;(s)  exp  (/^  X(n)  du)ds. 

0 s 

The  lemma  implies  that 

|E^(t)|  ^ |x^(t)|  + J^l A (s)  I I X^ (s) I exp  (/^  |A(u)|du)ds 


^ II  X^  II  {1  + /''|a(s)|  exp  (/^  I A (u)  I du)  ds } 


X.  II  • {1  + exp  ( 1 A (u)  I du)  ( /^-  ^ exp  (-/^  | A (u)  | du)  ds  )} 

^ 0 0 0 

X.||  • {1  + exp(/^|A (u)  |du)  (1  - exp (-/^l A (u)  I du) ) } 


= II  X . II  • exp  ( I A (u)  I du)  _<  II  X.  II  • exp  ( II  A II  t) 

II  II  20  E 

The  above  estimate  is  very  pessimistic  and  for  ||  A||  = 20,  e 
so  the  error  bound  will  be  very  large.  However,  we  can  use  this 
estimate  on  small  subintervals  as  initial  bound  for  the  error  and 
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then  improve  that  bound,  that  is:  for  u ^ t < u+h 

|E.(t)|  < |e.(u)|  + |x.(t)  - X.(u)l  + /*' 1 A (s)  I • I E . (s)  I ds 

hence 

IIe^  (•)  II  1 (1e.  (u)  I + ||x^  (•)  - X.  (u)  11  ) . (14) 

where  the  norms  are  computed  on  the  interval  [u,u+h].  If  h is  small 
the  above  estimate  is  not  a gross  estimate. 


Improving  the  bounds  on  |j  E,  )|  . 

Since  E^(t)  satisfies  equation  (13)  » one  has,  by  the  theorem 
in  section  (1.2) 

E^(t)  = Y(t)  Y ^ (s) A (s) X^ (s) ds  + X^(t) 

= (V(t)  + E (t))  /^(W(s)  + E (s))A(s)  X (s)ds  + X (t) 

1 0 ^ ^ ^ 

hence 

|e  (t)l  ^ (|v(t)l  + 1e  (t)l)  • /^(lw(s)|  + |e  (s) I ) 1 A(s) 1 |x(s) Ids 

^ 0 ^ 

+ |x^(t)]  . 

Similarly 

|e  (t)|  ;<(|w(t)|  + |e  (t)  |)/^(lv(s)  I + |e  (s)  I ) |a(s)  1 |X  (s)  Ids 

0 ^ 

+ lX2(t)l  . 

Using  the  above  estimates  we  arrive  at  the  following  algorithm: 
Notation; 

Let  0 = X,  < x_  . . . < X , = 1 be  a division  of  the  interval 

1 2 n+1 

[0,1] . 


h.  = X,  , - X, 

1 1+1  1 


a . ^ sup 

^ X . < t<  X 

1 1+1 


r^  ^ sup 


x.<t<x 
1 1+1 


lA(t) 1 


|X,  (t)  - X(x. ) I 
' 1 1 ' 
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^ sup  |x^(t)  - X^(x. ) 


x,<t<x.  , 
1 1+1 


‘2 ' ■ ' 2 i' 

|v(t)  I 


V,  ^ sup 
^ X . < t<  X , , 

i i+1 

w . ^ sup  I W (t)  I 

^ X . < t<  X . , 

i 1+1 

z , > sup  I X,  (t)  I 

z^  ^ sup  I X^  (t)  I . 

We  want  to  compute  such  that 

11 

i . > sup  f E,  (t) 

1 — 1 
X . < t<  X . , 

1 1+1 

£.  ^ sup  I £2^^^ 

X . < t<  X , , 

1-  - i+1 


define  ~ ~ 

Then  by  estimate  (14)  if 

n 

I = (8..  , + r.)e  ^ ^ 
1 1-1  1 

1°  = (£°  + ^.)e  ^ ^ 

1 1+1  1 

then 


8.?  ^ sup  I E (t) 

^ x.<t<x_  ^ 

1 1+1 


^ sup  I E2 

X , < t<  X . 

1 1 


also  if 


it  = (v,  + 8^)  y (w , + 8,)z,  a.h. 
1 11  3 3131 


i-1 


+ (v.  + i^)  (w.  + 8*?)z.  a.h.  + r. 
1 1 1 1111  1 


and 
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then 


-1  '0 

. = (w.  + 1.)  y (v . + J..)z,  a.h.  + 
1 3-  1.  , 1 

T = 1 


1-  (w  .+£*?*)  (v . +!l?)  z . a.h.  + r. 
11  11111  1 


Ie  (t)  I and  > sup  Ig  (t)  I 

1 x^t*^x  1 i.  ■”  X ^ t ^ X 2 ' ' ^ 

i-  — i+1  i—  — i+1 

So  we  set  1.  = min (£? , X, and  £.  = min(£*?,£l')  . We  also  would  like 

1 1 i 1 1 ii 

to  bound  ^2^^^  because  we  need  to  compute  a bound  on 

/^|(Y  ^(s))'  |ds.  Since  E'(t)  = A(t)  E (t)  + X'(t)  it  follows  that 


E-ll  1 11  A| 


^2 II  ^ 


The  above  estimates  enable  one  to  compute  good  bounds  on 

II  ^11  /^l  ^(s))'  Ids.  First  one  computes  numerically  the  values  of 

0 -1 

Y and  Y on  a sequence  of  points.  Then  interpolate  these  points 
by  a local  polynomial.  This  local  interpolation  scheme  will  give  V 
and  W . Then  using  the  above  algorithm  one  can  compute  rigorous 

bounds  on  II  y||  and  /^|(Y  ^(s))'  Ids. 

0 

y 

1.3.3  Computing  a Bound  on  ||  R ^ ||. 

In  the  previous  section  we  showed  how  one  can  compute 
6^  ^ |v(l)  - Y(l) I . Let 


D = Bi  + B^Vd) 


then 


If 


|r-d|  1 iB^lS^  = 62 


62  <1  then  R 


-1 


exists  and 


a) 


R~^l  < 


1-  D *6, 
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b) 


< 


- 1-  D 6, 


In  general  it  is  impossible  to  compute  D 

-1 


-1 


exactly.  Only  a good 


approximation  to  D 
that  I I -CD  I 


6 2 <1  then 


can  be  computed.  Let  C be  a matrix  such 

cl 


|d“^i  < 


- 1-6. 


Remarks : 

Computing  D ^ by  Gauss  elimination  might  not  yield  a very 
good  approximation  to  D However,  if  we  compute  a matrix  C such 
that  |i-Cd|  < 1,  this  approximation  can  be  improved  by  Newton's 
method.  The  condition  |i-Cd|  <1  is  enough  to  guarantee  that  the 
iteration  will  converge.  Therefore,  one  can  compute  an  approximate 
inverse  accurate  to  machine  precision. 

1.3.4  Computing  a bound  on  h = || 

An  obvious  bound  on  n is  ||  F '(Xg)  ^||"  II  » that  is,  the 

norm  of  the  Green's  function  times  the  residual.  If  the  norm  of  the 

Green's  function  is  not  very  large  the  above  estimate  is  good  enough. 

However,  this  is  an  overestimate  and  n is  usually  much  smaller. 

One  way  of  obtaining  better  bounds  on  n is  to  use  the  explicit 
form  of  the  Green's  function  in  order  to  compute  x^(t),  that  is; 

Let  v(t)  = Xj^(t)  - Xg(t),  r{t)  = xjt)-/^  f(s,Xg(s))ds  - Xq(0) 

then  by  the  theorem  in  section  (1.2) 

v(t)  = y(t)  /^(Y~^(s))^  r(s)ds  + Y(t)a  + r{t) 

0 

where  a is  given  by  equation  (10) . Therefore, 
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v(t)  = (V(t)  + E ft) ) { I ' (W  '(s)+E '(s) )r(s)ds  + a + E } + r(t) 
1^2  a 

0 

= V(t)  { W '(s) r(s)ds  + a}  + r(t) 

0 

+ E (t)  { /^  Y ^ (s) ' r (s)  ds  + a } 

^ 0 

+ V(t)  { E'(s)r(s)ds  + E } 

0 2 

where  a is  the  computed  value  of  a and  E^  = a-a.  Therefore 

II  '^11  II  V(")  ^ W ' (s)  r (s)  ds+a  II 

0 (15) 

+ II  E II  (/^|Y'^s)'|ds-  ||r  II  +11  dll)  + 1^/11  (IIe^II  -11  r||+  |e  1 + II  rll 
■^0 

In  orde^  to  compute  a,  one  can  use  the  following: 

Let  C be  the  computed  approximation  to  R ^ and  E^  = R ^-C. 
By  equation  (10)  (assuming  p=0) 

a = -R'^B^Yd)  y"^(s) ’r(s)ds  - R’^B^rd)  = 

= -(r“^B  -I)  y"^(s)' r(s)ds  - r"^B  rd)  = 

= -{[(C+E  )B  -I]/^(W'{s)+E '(s)  )r(s)ds  + (C+E  )B„r(l)} 
c 1 J 2 c 2 

= -{(CB  -I)/^  W'(s)r(s)ds  + C B rd)} 

^ 0 

- {e  B f W'(s)r(s)ds  + E'(s)r(s)ds]  + E B rd)} 

c 1 ^ 0 2 c 2 

= a + E 

a 

A bound  on  |E^j  is  given  by: 

|e  I 1 1e  1 • |bJ  (|  /^w  (s)r(s)ds|  +.  II  E'  11  • II  r|l  ) 

a c ± p 

+ Ie^I  • Ib^  rd)  I (16) 

Note  that  since  V,W,r  are  all  piecewise  polynomial  functions, 
V is  a piecewise  polynomial  function  and  all  the  necessary 
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j 


computations  can  be  carried  out.  The  ways  one  can  compute  bounds  on 

11.  IIe^  II  and  /^ly^(s)^lds  are  already  described  in  the 

c 1 Z Q 

previous  sections.  Therefore,  by  using  estimates  (15)  and  (16)  one 
can  compute  a bound  on  n . 

Another  possibility  for  computing  n is  to  solve  numerically 
the  linear  TPBVP: 
v'  = Av+r 

B^v(O)  + B^vd)  = 0 

where  r(t)  = f(t,XQ(t))  - x^(t)  and  A(t)  = f^(t,XQ(t)).  A bound 

on  II  ^]^”^qII  will  be  h = ||v||  + ||e'  (x^)  ^ ||  ’llqH  where 

q(t)  = v(0)  + A(s)v(s)ds  - v(t)  . 

0 

In  order  to  get  a better  bound  than  n = ||  F'  (x^)  ||  ||f(Xq)  || 

one  has  to  take  the  discretization  parameter  smaller  than  the  one 
used  in  the  original  equation  or  use  higher  order  method  (or  both) . 
Then,  hopefully,  ||q||  will  be  much  smaller  than  ||  F(Xq)||  . Since 
we  have  a IxDund  on  ||  F'  (x^)  ^ ||  we  know  how  small  ||  q ||  has  to 
be  in  order  to  get  a realistic  bound  on  n . 
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1 . 4 Bounding  the  Norm  of  Piecewise  Polynomial  Functions. 

In  order  to  compute  bounds  on  the  quantities  described  in  the 
proceeding  sections  one  needs  procedures  for  bounding  the  norm  of 
piecewise  polynomial  functions  on  a finite  interval.  The  polynomials 
may  have  interval  coefficients. 

In  this  section  we  are  considering  the  above  problem.  The  dis- 
cussion is  somewhat  elementary.  But  finding  effective  procedures 
for  bounding  the  range  of  piecewise  polynomial  functions  is  important 
enough  to  warrant  our  consideration.  We  wish  to  find  a good  balance 
between  speed  and  accuracy.  We  are  only  considering  robust  algo- 
rithms, that  is:  algorithms  that  are  guaranteed  not  to  give  erro- 


neous answers. 

We  begin  by  showing  that  bounding  the  range  of  interval  poly- 
nomials can  be  accomplished  by  bounding  the  range  of  at  most  4 poly- 


nomials. 

Let  P(x)  = 
val  coefficients. 


y A.x^  be  a polynomial  of  degree  n with  inter- 


We  wish  to  compute  ||  p||  on  the  interval  [a,b]  . 


Case  i)  0 a . 

In  that  case  let  q(x)  and  q(x)  be  defined  by 

n 

q{x)  = y sup  (A.  ) X 
i=0  ^ 

n 

q(x)  = y inf(A^)x^  . 

~ i=0  ^ _ 

Clearly  Vxc  [a,b)  g^Cx)  ^ P(x)  q{x)  and  q,^  are  the  best 
possible  functions  satisfying  the  above  inequality. 
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Case  ii) 


b < 0 


Define 


n 


q(x)  = 1 sup  ((-l)^A.)x 
i=0  ^ 


n 


1 1 
X 


£(x)  = J inf  ( (~1)  A ) 
i=0 

again.  x £ [a,b]  £(x)  ;^P(x)  ^ q(x) 

Case  iii)  a < 0 < b 

Define 


i(x)  = < 


a(x) 


n 

sup  ( (-1) ^A  , ) X 
1 

I 

i=0 

H- 

II 

O 

sup (A^) x^ 

n 

inf ((-l)^A,)x^ 
z 

I 

i=0 

n 

I 

i=0 

inf  (A . ) 

1 

a < X < 0 


0 < X < b 


a < X < 0 


0 < X < b 


Therefore,  by  bounding, at  most,  4 polynomials,  one  can  compute  bounds 
on  the  norm  of  P(x).  In  our  programs  all  polynomials  are  defined 
on  the  interval  [0,h]  where  h is  the  size  of  the  subinterval  . 
Therefore,  only  case  i)  applies.  In  practice  we  resort  to  the  above 
reduction  only  if  the  width  of  the  coefficients  exceed  a prespecified 
tolerance. 


1.4.1  Computing  a "Rough"  Bound  on  Piecewise  Polynomial  Function. 

If  the  interval  on  which  the  polynomial  pieces  are  defined  are 
small  and  the  piecewise  polynomial  function  itself  is  not  very  small, 
one  can  use  the  following  effective  procedure: 

Assume  that  the  polynomial  p(x)  is  defined  on  the  interval 

[0,h]  . 


I 

b 
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a)  Determine  upper  and  lower  bounds  on  P'(x); 


q ^ P'  (x)  ^2 


X e [0,h] . 

One  way  to  compute  q,2  is  as  follows: 
n-1 

if  P' (x)  = y b.x^  let  b.  = max(0,b.)  b.  = min(0,b.) 

^ T T ”11  1 


n-1 

n-1 

i = b^  + 1 b.h^  , 

a = 

b„  + 1 b.h^ 

0 -1 

i=l 

1=1 

b)  if  P*  (x)  > 0 or 

P'  (x) 

< 0 determine  ||  p|| 

by  evaluating 

P(x)  at  the  appropriate  end  point. 


c)  If  2”^0<q,  then  set 

II  p||  = max(  |p(0)+hq|  , |p(0)+h£|  ) . 

A better  bound  can  be  computed  by  making  use  of  the  value  of  P(h). 


Set  II  P II  = max  ( I P (0) +Xhq|  , |p(0)+ (1-X ) h^l  ) where 

P(h)  - P(0)  - hg 

^ — - — • 

hq  - hg 

It  is  our  experience  that  the  above  procedure  gives  bounds  with 
enough  accuracy  for  bounding  piecewise  polynomial  functions  which  are 
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not  small  and  the  right  order  of  accuracy  of  the  bounds  on  the 

norms  of  piecewise  polynomial  functions  that  are  very  small. 

We  find  the  above  procedure  to  be  adequate  for  computing  bounds 

on  II  a||  , II  v||,  II  w||  and  for  computing  bounds  on  Hw-I-J  WAds|| 

0 

and  ||v-I-/AVds  ||  . The  reason  being  that  we  only  need  to  compute 
0 

bounds  that  agree  with  the  best  possible  bounds  up  to  one  or  two 

digits  and  we  do  not  need  the  best  possible  bounds.  In  many  cases 

the  above  procedure  is  adequate  for  bounding 

II  v(t)  - v(0)  - f (x,v(s)  )ds||  . 

0 

1.4.2  Computing  a Good  Bound  on  the  Norm. 

In  order  to  compute  a good  bound  on  the  norm  of  a polynomial 

P(x) , one  can  bracket  the  real  zeros  of  P' (x)  inside  very  small 

intervals.  Then  one  uses  step  c)  in  the  previous  method  in  order 

to  bound  the  norm.  One  well-known  method  for  bracketing  the  real 

zeros  of  a piolynomial  is  the  Sturm  Sequences  method  (See  for  excimple 

Isaccson  and  Keller  [ 19] ) . Another  way  to  bracket  is  to  use  an 

algorithm  suggested  by  Dussel  and  Schmitt  [13].  First, one  finds  an 

(r) 

integer  r such  that  the  rth  derivative  of  P,  P does  not  have 

(r—l ) 

a zero  in  the  interval.  Using  the  fact  that  P (x)  can  have  at 

(r-1) 

most  1 zero  one  can  easily  find  out  if  P does  have  a zero  or 

(k-1) 

not.  If  it  does,  the  zero  can  easily  be  bracketed  since  P is 

( ir“lc) 

monotone.  In  the  kth  step,  assuming  that  the  zeros  of  P (x) 

{ 3r”K) 

are  bracketed,  one  can  use  the  sign  of  P (x)  on  each  subinter- 

( ir“}c“l ) 

val  in  order  to  bracket  the  zeros  of  P (x) . 

We  have  chosen  to  follow  the  idea  of  Dussel  and  Schmitt,  but 
instead  of  using  an  interval  version  of  the  Bisection  Method  as  they 
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do,  we  use  an  interval  version  of  Modified  Regula  Falsi  algorithm. 
(For  details  of  Modified  Regula  Falsi  algorithm  see  Conte  and 
de  Boor  [9  ]).  We  use  the  above  method  only  on  intervals  of  which 
the  function  is  tionotone.  In  the  case  of  multiple  zeros  we  are  using 
a simple  gross  bound.  In  that  case  the  function  is  very  flat  and  a 
good  bound  on  the  zeros  of  the  derivative  is  not  needed.  A careful 
examination  of  all  possible  cases  shows  that  the  cibove  procedure 
gives  guaranteed  lower  and  upper  bounds  for  the  zeros.  The  details 
of  the  above  algorithm  are  given  in  the  next  section. 

One  last  consideration.  The  intervals  on  which  each  polynomial 
is  defined  are  small.  Therefore,  the  contribution  to  the  norm  from 
the  high  powers  is  very  small.  One  can  bound  the  contribution  from 
the  low  powers  and  then  add  the  contribution  from  the  absolute  value 
of  the  high  powers.  This  way  the  degree  of  the  polynomials  that  one 
has  to  bound  is  reduced.  Note  that  it  is  not  essential  to  compute 
the  bounds  very  accurately.  Accuracy  of  two- three  digits  is  good 
enough. 

1.4.3  On  Bounding  the  Range  of  a Polynomial. 

k 

Let  I = [a,bl  , a 0,  be  a finite  interval  and  P(x)  = a.x 

i=0  ^ 

a kth  degree  polynomial.  In  order  to  bound  the  range  of  P on  I 

we  bracket  the  zeros  of  P' , that  is:  We  find  a sequence  of  points 

a<x,  <x.,...<x  = b such  that  on  each  interval  [x.  ,x.  ,] 

— 1 2 n 1 i+i 

P' (x)  > 0,  P' (x)  <0  or  P' (x)  may  have  a zero.  We  try  to  make 
the  interval  on  which  P*  (x)  may  have  a zero  as  small  as  possible. 
Once  we  have  bracketed  the  zeros  of  P'  inside  small  intervals, 
we  use  the  procedure  described  in  the  preceding  section  in  order  to 
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bound  the  norm  of  the  polynomial. 

In  order  to  bracket  the  zeros  of  P'  we  use  recursively  the 
following  elementary  fact  (Rail's  theorem). 

fact : If  f (x)  has  no  zeros  in  an  interval,  then  f has  at  most 

one  zero  in  that  interval. 

We  start  by  finding  an  integer  j as  small  as  possible  such 
that  (x)  has  no  zeros  in  I . Clea  ly  there  is  at  least  one 


such  j namely  j = k. 

If  j = 1,  then  we  are  done.  That  is  P' (x)  has  no  zeros  in 

I . Otherwise,  we  assume  we  have  a sequence  of  points  a = x^  < x^  ^ 

...  < X = b and  a list  of  integers  S,,S.,,  ...,  S . such  that 
n 1 2 n-1 

S.  € {-1,0,1}  and 
1 


If 

s. 

= 1, 

p^3) 

(x) 

> 0 

on 

[x, ,X.  - ] . 

1 

1 1+1 

If 

s. 

= -1 

p(3) 

(x) 

< 0 

on 

(x. ,X,  - ] . 

1 

t 

1 1+1 

If 

s . 

1 

= 0, 

p(j) 

(X) 

may 

have 

a zero  in 

^i"'i+l'  • 

We  always  try  to  make  the  intervals  with  possible  zero  smaller  than 
some  tolerance  e . 

Given  such  a sequence  and  a list  of  integers,  we  compute  a new 
sequence  of  points  a = y^  < y^  < • • • ~ ^ list  of  inte- 


m-i 

such 

that 

V.  = 1 

1 

if 

p(j-l) 

(X)  > 0 

on 

^^I'^i+i’ 

V.  = -1 
1 

if 

p(3-l) 

(X)  < 0 

on 

f^i'^i+i’ 

o 

It 

•H 

> 

if 

p(D-l) 

may  have 

a zero  in  [y^  * 

Clearly 

at  most 

k-1 

such  steps 

are 

needed  in  order  to 

zeros  of  P' (x)  on  I 
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We  now  describe  the  procedure  by  which  we  construct  the  se- 
quences y,  : V with  the  aid  of  the  se- 

12  m 12  m-l 

quences  x, ; S S . .,S  For  convenience  we  denote 

12  n 1 2 n-1 

by  f and  p ^ by  f . 

We  scan  the  interval  I from  left  to  right.  On  each  subinter- 
val [x.,x.  ,]  we  determine  if  f has  a zero.  If  it  does  not  we 
1 i-tl 

proceed  to  the  next  subinterval.  If  it  does  have  a possible  zero, 
then  we  bracket  that  zero  in  an  interval  as  small  as  possible  (up  to 
a tolerance  e) . In  order  to  take  round-off  error  into  account,  f 
is  evaluated  in  Interval  Arithmetic.  So  at  any  point  x that  we 
evaluate  f we  get  two  numbers  f(x),  £(x)  such  that 
f (x)  > f (X)  ^ £(x)  . 

One  has  to  consider  the  following  different  cases: 


Case 

_i 

s.=l 

1 

• 

In  that  case 

f is  strictly  increasing. 

Therefore,  f 

has  a 

zero 

in 

[x^  ,x 

i-H^ 

if  and  only  if  f(x.)  < 0 
1 — 

and  f (x.  , - ) 
1-1-1 

1 0 . 

Case 

1 . a 

If 

f(x.) 
— 1 

> 0 

or  f (x.  , ) < 0 f does 
i-i-l 

not  have  a zero  in 

[X.,: 

’'i+i’ 

• 

Case 

l.b 

If 

f (x.  ) 
1 

< 0 

and  ^ ® ' then 

f has  a zero 

in  the 

open 

inte rval 

(x. 

1 

,Xi_^l)  . 

We  can  trap  that  zero  inside  a small  interval  by  a disection 
method  (Modified  Regula  Falsi)  which  we  will  describe  later. 
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Case  1 . c 


f(x.)  < 0 and  f(x.  ,)  >0  but  f(x.)  > 0 
— 1 — — 1+1  1 — 


There  are  two  different  possibilities. 


1)  i = 1 in  that  case  set  y 


V = 0.  Since  f(x.)  < 0 
1 — 1 — 


and  f(x,  -)  > 0 we  can  use  the  disection  method  to  find  y such 
that  jf(y)  > 0 and  y is  as  small  as  possible  (that  is:  y is  a 
machine  representable  number  and  f^(y-e)  ^ 0).  Set  = 'i,  ^ 

and  go  to  the  next  interval  since  f has  no  zeros  in  the  interval 

[y,x.+i]  . 

2)  i > 1 In  that  case  we  have  y^  such  that  f might  have  a 
zero  to  the  right  of  y^  and  = 0.  Again  we  can  find  y as 

small  as  possible  by  the  disection  method  such  that  f(y)  > 0 . Set 
y^^j^  = y , ^ 9°  next  interval.  Note  that  if 

f(x. ) =0  we  can  take  y = x^  , where  xT  is  the  smallest  machine 
representable  number  greater  than  x^  because  ^(x^)  ^ f (x J and 
f is  strictly  increasing  on 
Case  1 ■ d 

f(x.)<  0 and  f(x..,)  > 0 but  f(x,  ,)  < 0 
1 1+1  — — 1+1  — 
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1)  If  I = 1 set  ~ find  y as  large  as  pos- 
sible such  that  f(y)  < 0 . Set  y^  = y and  ^^2  ~ ' 

2)  If  i > 1 we  have  y.  such  that  f(y.)  < 0 and  V,  = -1.  Find 


0. 


j+1 

Case  l.e 


f(x,)  < 0 < f(x,)  and  f(x.  ) < 0 < f{x.  ). 

— 1 — — 1 — 1+1  — — 1+1 

If  i = 1 then  set  y^^  = = 0 and  go  to  the  next  interval. 

If  i > 0 then  go  to  the  next  interval.  (Remark:  this  is  a very 

unlikely  possibility) . 

Case  2 S.  = -1  . 

1 

This  case  is  almost  the  same  as  Case  1.  The  only  difference  is 

that  the  function  is  strictly  decreasing. 

Case  3 S.  = 0. 

1 

Case  3. a 


There  is  no  possible  zero  in  the  previous  interval  or  i = 1. 

In  that  case  we  compute  f([x^,x^^^])  by  interval  arithmetic. 

If  there  is  no  zero  we  continue  with  the  next  interval.  If  there  is 

a zero  in  f([x. ,x.  ,])  then 
1 1+1 


If 

i = 

1 then  y.  = x,  , 

o 

11 

and  if  0 4 [^(x^),  fCx^)] 

^2 

= ’‘2 

, = sign  (f  (x^) ) . 

otherwise  go  to  the  next  interval. 

If 

i > 

0 set  y , , = X,  , 

V 

0 . 

■']+l  1 

3+1 

If 

0 4 

[f(x.^^f(x.^^]  set 

= Vl'  ''j+2  = 

else  go  to  the  next  interval. 
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Case  3 . b 


There  is  a possible  zero  in  the  previous  interval. 

If  0 4 [f(x  ),f(x.  )]  set  y = x and  V,  = 

— 1+1  1+1  j+1  1+1  3+1 

sign (f (x^^^) ) , otherwise  go  to  the  next  interval. 

Remark : The  above  procedure  does  not  give  the  tightest  possible 

intervals  in  cases  3. a and  3.b.  However,  since  in  these  cases  the 
function  is  very  flat  in  that  neighborhood,  one  does  not  really  need 
a very  tight  bound  on  the  zeros  of  the  derivative.  In  [13]  Dussel 
and  Schmitt  describe  how  one  can  use  an  interval  version  of  the  bi- 
section method  to  get  tight  bounds  on  the  zeros  even  if  the  function 
is  not  monotone.  One  can  use  that  procedure  if  one  wishes  to  . 

1.4.4  The  Dissection  Algorithm. 

The  algorithm  we  use  is  a Modified  Regula  Falsi  algorithm  (see 
Conte  and  de  Boor  [ 9 ] ) . We  use  two  versions : one  for  functions  that 
take  interval  values,  and  the  other  for  functions  that  take  real 
values.  Both  algorithms  differ  from  each  other  by  minor  details. 

We  assume  that  we  are  given  an  interval  [a,b]  such  that 
f(a)  • t(b)  < 0 . At  each  step  we  have  FA  and  FB  such  that 
FA'FB  < 0.  We  compute  X by:  X = 1/ (1-FA/FB) . Clearly  0 f.  1 ^ 1 • 

In  order  to  avoid  problems  of  convergence  and  overflow-underflow, 
if  X < 6 we  set  X = 6 or  if  1-X  <6  we  set  X = 1-6.  6 is  a 

given  small  number,  6 > 0.  After  we  determine  0 < X <1/  we  conpute 
a point  X in  the  interval  (a,b)  by  x = X*a  + (1-X) • b.  If 
0 e f(x),  then  we  have  found  a zero,  otherwise  depending  on  the  sign 
of  f (x) , we  replace  a or  b by  x and  FA  or  FB  by  f(x). 
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If  a (or  b)  had  been  replaced  in  the  previous  iteration  and  is 
being  replaced  again  we  set  FB  = FB/2  (or  FA  = FA/2) . The  above 
step  prevents  the  possibility  of  approaching  the  zero  only  from  one 
side.  If  b-a  ^ z then  we  are  done.  Otherwise,  we  confute  new  A 

and  so  on. ..  .In  the  case  where  f is  an  interval  valued  func- 

tion, If  0 e f(x),  then  we  check  if  0 4 f (x  - ^)  and 
0 4 f (x  + j)  ; if  this  is  true  we  take  (x  - j,  x + j)  to  be  the 
desired  interval. 

If  0 e f(x  - ^)  or  0 e f(x  + ^)  then  we  use  the  same  proce- 
dure for  f and  ^ in  order  to  find  the  intervals  containing  zero 
of  f and  ^ . Using  these  intervals  and  the  sign  of  f we  deter- 
mine an  interval  containing  the  zero  of  f . 

In  case  f(x)  =0,  [x  ,x^]  is  an  interval  containing  a zero 

of  f and  the  same  is  true  for  f . (x  is  the  biggest  machine 

number  smaller  than  x and  x^  is  the  smallest  machine  number 

greater  than  x).  The  above  algorithm  is  linearly  convergent  since 

the  interval  is  reduced  at  least  by  a factor  of  1-6  at  each 

iteration  and  at  most  by  a factor  of  6 . However,  if  6 is  small 
-2  -3 

(say  6 = 10  or  10  ) the  algorithm  behaves  like  a quadratic- 

ally  convergent  algorithm.  In  some  sense  this  algorithm  is  in 
between  the  Bisection  method  (6  = j)  and  the  modified  Regula  Falsi 
(6=0). 
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CHAPTER  2 


A POSTERIORI  ERROR  BOUNDS  - THE  IMPLEMENTATION 
2.1  On  *-.he  Numerical  Procedures. 

In  this  section  we  describe  the  nxjmerical  procedures  we  use  in 
order  to  compute  upper  bounds  on  the  constants  of  Kantorovich's 
theorem. 

Let  us  assume  that  we  are  given  a numerical  solution  to  equa- 
tion (1)  computed  at  a sequence  of  points  a = x,  < < . . . < x , 

1 2 n+i 

= b.  As  an  initial  guess  we  take  the  piecewise  polynomial  function 
defined  as  follows:  On  each  subinterval  ^ ~ l,...,n-l 

x^Ct)  is  the  sixth  order  polynomial  that  interpolates  the  solution 
and  its  derivative  at  the  points  x.,  x.  ,,x.  The  values  of  the 

derivatives  are  computed  via  the  differential  equation.  On  the  in- 
terval ^^n'^n+1^'  ^0^^^  sixth  order  polynomial  interpolat- 

ing the  solution  and  it's  derivative  at  x ,,  x , x 

n-1  n n+1 

First  the  residual  is  computed,  that  is: 

r(t)  = x^Ct)  - x^(a)  - f (s  ,x^  (s)  )ds . 

a 

On  each  subinterval  ^ polynomial.  Since  f(s,y)  is  a 

polynomial  in  s and  y , on  each  of  the  subintervals  r(t)  is 
also  a polynomial.  Using  the  polynomial  coefficients  of  r(t)  and 
the  algorithm  described  in  section  (1.4.2)  a bound  on  |1  r||  is 
computed . 

Second,  the  fundamental  and  inverse  fundamental  solutions  are  nu- 
merically computed.  One  has  to  solve  2n  Initial  Value  Problems: 

-1  T 

n for  y(t)  and  n for  (Y  (t) ) . The  values  of  Y(t)  and 
-1  T 

(Y  (t))  are  computed  at  the  same  sequence  of  points  that 


was  computed. 

We  take  V(t)  to  be  the  piecewise  polynomial  function  approxi- 
mating Y(t),  and  W(t)  to  be  the  piecewise  polynomial  function 
approximating  \ ^(t).  V(t)  and  W(t)  are  computed  with  the  same 
interpolation  scheme  used  to  construct  XQ(t). 

The  algorithm  described  in  section  (1.3.2)  is  used  to  confute 

bounds  on  11  E,  ||  # II  E , 1|  Y |1  , / | Y ^ (s)  ' | ds  and  1 E (1)  | . 

a 

The  bounds  on  lix^||  , 11x2  H (of  equation  13a,  13b)  as  well  as 
the  bounds  on  |1  vl|  , ||  W 1|,  H W'  1|  and  1|  a|1  are  computed  using 
the  rough  bound  described  in  section  (1.3.2)  . 

Third,  R = + B2  V (b)  and  C = R ^ are  computed.  Using 

the  bounds  on  |e(1)1  and  Icl  bounds  on  1r  ^1,  1r  ^-cl 
max ( 1r~^B  1 , |r”^B  -I  I ) and  1r  ^B  1 are  computed  (see  section 

X -L  ^ 

(1.3.3))  . 


The  facts  that  the  problem  has  separated  boundary  conditions 
and  that 


r I 0 -i  „ f 0 0 

®1  ~ ^ 0 0 ^ ' ®2  * I 0 ^ 


implies  that: 

a) 

IbJ 

= IB2I 

= 1 

b) 

if 

Y(b)  = 

f a 3 1 
'■  r 6 *' 

where  a,  6,  y,  <5 

are 

n/2  X 

n/2  matrices,  then 
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Therefore,  max(lR  ^B^l,  |r  = [S  ^u|  + 1 and  jR  1 = 

|S  ^1  . Using  the  above  bounds  the  bound  on  ||  F'  (x^)  ^ ||  is  com- 

puted by  formula  (12). 

The  bound  on  ||  h r given  by 

n = II  yII  ( |y  ^(s)  Ids  - ( |B  ^a|  + 1)  + I B ^ I ) II  r II  + ||  r||. 
a 

Since  the  problem  has  simple  boundary  conditions  one  makes  sure  the 
initial  guess  satisfies  the  boundary  conditions.  Therefore, 

one  can  assume  p = 0.  The  above  bound  for  n is  smaller  than  the 
bound  (I  F'  (Xq)  ^ i|  • II  rjl  . 

2 . 2 Function  Representations . 

The  numerical  solution  to  equation  (1),  the  fundamental  solu- 
tion and  the  inverse  fundamental  solution  are  continuous  functions. 
However,  these  functions  have  a discrete  representation  inside  the 
computer.  We  are  using  two  types  of  representations: 

a)  Point-value  representation:  In  this  representation  the 

function  is  represented  by  a sequence  of  points  (the  knot  sequence) 

a = x,  <x^<...<x  = b and  a sequence  of  values  y,  ,..-,y 

12  n I n 

whicli  are  the  double  precision  values  of  the  function  at  x, ,...,x  . 

In 

We  use  the  same  knot  sequence  to  represent  x^{t),  the  approximate 
solution  of  equation  (1),  V{t),  the  approximate  fundamental  solu- 

tion and  W(t),  the  approximate  inverse  fundamental  solution. 

b)  Piecewise  polynomial  representation:  In  this  representa- 

tion the  function  is  represented  by  the  knot  sequence  and  a 6xn 
array  of  numbers.  On  each  subinterval  function  is 

a 6th  order  polynomial  p(t)=y  a.  . (t-x.)^.  This  polynomial  inter- 

j=0  3'"  " 

polates  the  values  of  the  function  and  its  derivative  at  the  points 
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Since  computers  compute  with  a finite,  discrete  set  of  numbers, 
it  is  impossible  in  general,  to  compute  the  exact  coefficients  of 
the  polynomials.  Moreover,  the  coefficients,  usually,  are  not  ma- 
chine representable  numbers.  On  the  other  hand  the  initial  guess 
XpCt),  V (t)  and  W(t)  are  to  be  continuous  functions.  However,  one 
cannot  guarantee  that  the  computed  functions  are  continuous.  To 
overcome  this  diffuculty  we  are  using  interval  arithmetic  (see  [30]). 
Instead  of  computing  with  one  floating  point  number  the  computation 
is  done  with  two.  One  of  the  numbers  is  guaranteed  to  be  smaller 
or  equal  to  the  number  we  need  and  the  other,  greater  or  equal  to 
it.  Therefore, one  actually  has  two  sets  of  polynomial  coefficients. 
The  first  is  a machine  representable  set  of  numbers.  Each  of  them 
is  guaranteed  to  be  greater  or  equal  to  the  corresponding  exact  co- 
efficient. Similarly  each  floating  point  number  in  the  other  set 
is  guaranteed  to  be  smaller  or  equal  to  the  corresponding  exact  co- 
efficient. Therefore,  one  has  two  polynomia?s;  one  which  at  any 
point  is  greater  or  equal  to  the  exact  polynomial,  and  the  other  is 
smaller  or  equal  to  it.  Although  it  is  impossible  to  compute  the 
exact  polynomial  itself,  by  bounding  the  range  of  these  two  poly- 
nomials one  gets  a bound  on  the  range  of  the  exact  polynomial . 


2 .3  The  Interval  Arithmetic  package. 

As  was  said  before  we  use  interval  arithmetic  in  order  to  take 
the  effe-t  of  round-off  error  into  acco\int. 

The  interval  arithmetic  package  we  have  constructed  is  based 
on  the  Multiple  Precision  Package  (See  Crary  [lo] ) . The  interval 
precision  package  was  designed  to  be  used  with  the  AUGMENT  precom- 
piler. All  subroutines  that  use  Interval  or  Multiple  arithmetic  are 
written  in  "extended"  FORTRAN  (See  [11]),  and  are  translated  by 
AUGMENT  to  standard  FORTRAN.  For  the  use  of  AUGMENT  see  Crary  [11]. 
The  package  is  similar  to  the  one  written  by  Yohe  [44].  However, 
the  package  is  a "scaled  down"  version;  only  routines  that  are 
needed  were  constructed.  Yohe's  package  was  not  used  because  we 
needed  more  precision  than  that  package  provides.  Yohe's  package 
is  written  in  single  precision  arithmetic  (approximately  8 decimal 
digits)  . Our  package  uses  4 words  per  floating  point  numlser. 

(About  32  decimal  digits.)  Therefore,  each  interval  number  takes  8 
words . 

The  subroutines  and  functions  available  are: 

1)  The  arithmetic  operations:  +,  -,*,/. 

2)  Comparison  operators:  .LT.  , .LE. , .GT.,  .GE. , .EQ.,  .NE. . 

3)  Functions:  INF,  SUP,  ABS,  TSIGN,  COMPOS,  MAX.  INF,  SUP, 
ABS  and  COMPOS  are  the  usual  ones  TSIGN  is  defined  by: 

p if  0<a 

TSIGN(  [a,b]  )=/o  if  a_<0£b 

1-1  if  b < 0 
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and  MAX  is  defined  by: 


MAX ( La ,b] , [c,d] ) = [f,f]  where  f = max(b,d). 

4)  Conversion  functions: 

There  are  conversion  functions  from  Integer,  Double  Precision, 
and  Multiple  to  Interval  and  from  Interval  to  Double  percision  and 
Integer.  Type  conversion  has  to  be  explicitly  invoked  and  is  not 
done  automatically.  That  is,  the  statement  AI  = D where  AI  is 
of  type  Interval  and  D is  Double  precision , is  not  legal.  The 
function  CTX ( • ) converts  its  argument  to  Interval,  CXTD  converts 
Interval  to  double  precision  and  CXTI  converts  Interval  to  Inte- 
ger. Since  the  Multiple  precision  package  is  used  to  construct  the 
interval  package  the  description  deck  of  the  multiple  package  has  to 
be  submitted  to  AUGMENT  in  front  of  the  Interval  description  deck. 
The  description  decks  of  Multiple  and  Interval  are  given  in  the 
Appendix . 
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2 . 4 The  Subroutine  Package . 


2.4.4  General  Information. 

This  package  of  subroutines  enables  one  to  bound  ||  x^-x^H 
and  II  F'  ^||.  Most  of  the  computation  is  done  in  multiple  and 

interval  arithmetic.  However,  the  user  does  not  have  to  be  very 
knowledgeable  about  the  inner  working  of  the  multiple  and  interval 
packages.  The  main  subroutine  APOSTR  does  not  have  multiple  or 
interval  arguments  in  its  calling  sequence. 

There  are  five,  user  supplied,  problem  dependent  subroutines. 
Two  subroutines  out  of  the  five  use  Interval  arithmetic.  In  order 
to  be  able  to  write  these  routines  the  user  has  to  know  how  to  use 
AUGMENT  (see  [11])  and  how  to  use  a few  polynomial  manipulation  sub- 
routines. The  polynomial  subroutines,  the  interval  package  and  the 
description  decks  are  provided  with  the  package.  The  underlying 
arithmetic  is  provided  by  the  multiple  precision  package  and  is 
available  on  MRC*LIB. , that  is  Mathematics  Research  Center's  re- 
locatable library.  Since  the  multiple  package  is  written  in  UNIVAC 
1100  series  assembler,  the  package  is  not  portable  and  no  attempt 
was  made  to  make  it  portable.  We  tried  to  make  the  package  applic- 
able to  a large  set  of  problems.  However,  the  package  was  designed 
as  an  experimental  tool  to  test  some  of  our  ideas.  It  is  not  and 
was  not  meant  to  be  a general  purpose  production  code. 

We  now  give  a short  description  of  the  subroutines  in  the  pack- 
age. A complete  listing  of  the  subroutines,  the  extended  FORTRAN 
source  code  as  well  as  the  Standard  FORTRAN  code  produced  by 
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4. 


AUGMENT  is  given  as  an  appendix  on  a microfiche  card. 

2.4.2  Problem  Independent  Subroutines . 

APOSTR : This  is  the  main  subroutine.  It  sets  up  the  work  space 

for  each  of  the  subroutines,  calls  the  subroutines  that  perform  the 
different  parts  of  the  algorithm  and  computes  the  final  bounds. 

See  calling  sequence  description  and  explanations  in  the  next  sec- 
tion. 1 

INLIZ:  Initialization  routine  sets  up  some  needed  constants. 

SETKNT : This  subroutine  converts  the  double  precision  kn-^t  se- 

quence to  interval  valued  knot  sequence. 

GREENF : This  subroutine  puts  the  values  of  the  solution  into  com- 

mon block  /SPLCFF/  and  also  computes  the  values  of  Y(t)  and 
-1  T 

(Y  (t) ) at  the  knot  sequence.  It  uses  subroutine  DDESUB  in 
order  to  solve  the  initial  value  problems. 

CMPORS : This  subroutine  computes  a bound  on  the  norm  of  the  resid- 

ual, that  is:  A bound  on  ||  ^c^Ct)  - /^f  (s  , x^  (s ) ds  - x^  (a)  |1  . 

a 

BDGRFC : This  subroutine  implements  the  algorithm  described  in 

section  (1.3.2).  It  computes  bounds  on  ||Y(t)|l  , /^  | Y ^ (s)  ' | ds,  ||  A ||  , 

a 

II  eJI  , II  E^li  and  lE^(b)  |. 

CMRINV ; This  subroutine  confutes  an  approximation  to  R ^ where 
R = B^  + B^YCb).  The  approximation  is  first  computed  by  Gauss 
elimination  and  then  improved  by  Newton's  method.  The  subroutine 
computes  bounds  on  ||  R ^ ||  , ||  B^^R  ^ 1|  and  ||  B^R  ^ |1  using  the 
method  described  in  section  (L3.3)(see  also  section  (21)). 

^Only  the  MRC  report  has  the  microfiche  card  for  obvious  reasons.  A 
copy  of  the  microfiche  card  can  be  obtained  from  the  author. 
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T.F.QRMR  • This  subroutine  computes  the  residual  of  V , that  is: 

V(t)  - /''A  (s)  V (s)  ds  - V(a).  (Recall  that  V(t)  is  the  piecewise 
a 

polynomial  approximation  to  Y(t).)  LEQRMR  is  used  by  subroutine 
BDGRFC . 

LQERMR : This  subroutine  computes  the  residual  of  W , that  is 

W(t)  - W(a)  + W(s)A(s)ds.  (Recall  that  W(t)  is  the  piecewise 
^ -1 

polynomial  approximation  to  Y (t) ) . LQERMR  is  used  by  BDGRFC. 
CBSTPP : This  subroutine  converts  point  value  representation  of  a 

function  to  its  piecewise  polynomial  representation.  It  confutes 
the  coefficients  of  the  polynomial  interpolating  the  function  and 
it's  derivative  at  3 consecutive  points. 

PPNRM:  This  subroutine  computes  a tight  bound  on  the  norm  of  a 
polynomial  in  a given  interval . It  uses  subroutine  PNORM  which 
is  the  implementation  of  the  algorithm  described  in  section  (1.4.2). 
PNORM:  This  subroutine  implements  the  algorithm  described  in  sec- 

tion (1.4.2).  It  isolates  the  zeros  of  the  polynomial's  derivative 
inside  very  small  intervals  and  then  bounds  the  norm  of  the  poly- 
nomial. It  uses  subroutine  LOCATE  in  order  to  bracket  the  zeros  of 
the  polynomial's  derivatives  and  subroutine  NORMl  in  order  to  com- 
pute bounds  on  the  norm  of  the  polynomial. 

NROMl : This  subroutine  is  given  a sequence  of  intervals  containing 

the  zeros  of  the  polynomial' s derivative  using  this  information  it 
bounds  the  polynomial's  norm. 

LOCATE ; This  subroutine  is  part  of  the  algorithm  for  getting  a tight 
bound  on  the  norm  of  a polynomial.  This  subroutine  gets  a sequence 
of  points  where  the  kth  derivative  of  the  polynomial  is  either 
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positive  or  negative  or  may  have  a zero.  Einploying  the  method 


described  in  section  (1.4.3)  subroutine  LOCATE  computes  a sequence 
of  points  where  the  (K-l)th  derivative  is  positive  or  negative  or 
may  have  a zero.  The  intervals  that  could  possibly  contain  a zero 
are  made  smaller  than  a given  tolerance. 

BISECT:  This  is  the  dissection  algorithm  described  in  section 

(1.4.4).  This  subroutine  brackets  the  zeros  of  sup  or  inf  of 
an  interval  valued  polynomial.  The  polynomial  is  assumed  to  have  a 
zero  in  the  given  interval.  It  is  also  assumed  to  be  monotone. 
BSCINT : The  same  as  BISECT  only  it  is  the  interval  version  of  the 
dissection  algorithm. 

PMXNRM : This  subroutine  computes  a bound  on  the  uniform  norm  of  a 

polynomial  matrix.  The  bound  is  given  by:  max  'l  |(  a.  . ||  . The 

1 li^n  j=l 

are  computed  by  the  "rough  bound" 
method  described  in  section  (1.4.1). 

a 

UBEXP : This  subroutine  computes  a bound  on  e where  a is  a 

positive  interval  numiser. 

VCTINT : This  subroutine  integrates  piecewise  polynomial  functions. 

X . 

That  is:  it  is  given  / u(s)ds  and  the  polynomial  coefficients 

a 

of  u(t),  x^  ^ t ^ ^i+1'  computes  the  polynomial  coeffi- 
cients of  u(s)ds,  X,  < t < X.  , . u is  a vector  valued  func- 

a 

tion . 

PPVLUl : This  subroutine  computes  the  double  precision  value  of  the 

solution  at  any  given  point. 

DIVDIF : Divided  differences  subroutine.  It  is  used  by  PPVLUl  to 

compute  the  interpolating  polynomials. 


bounds  on  the  terms  ||  a.  . 1| 
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IDVDIF: 


The  same  as  DIVDIF  only  in  interval  arithmetic.  It  is 


used  by  CBSTPP  in  order  to  compute  the  polynomial's  coefficients. 
DDESUB : This  subroutine  solves  initial  value  problems  ii;  double 

precision.  This  subroutine  uses  Gragg's  modified  midpoint  rule 
with  Bailirsch-Store ' s rational  extrapolation  method.  It  is  es- 
sentially the  same  subroutine  as  the  one  provided  by  Madison 
Academic  Computing  Center.  A modification  was  made  to  enable  the 
routine  to  store  the  solution's  values  at  any  given  sequence  of 
points  not  necessarily  equally  spaced.  For  details  about  the  use 
of  the  subroutine  see  [49] . 

2.4.3  Polynomial  Manipulation  Subroutines. 

In  order  to  facilitate  the  manipulation  of  polynomials  several 
subroutines  were  constructed.  All  polynomials  have  interval  coef- 
ficients. The  following  siJDroutines  are  provided: 

SUBROUTINE  PLYNEG (K, A,B) 

This  subroutine  computes  B = -A,  A and  B are  Kth  order 
polynomials . 

SUBROUTINE  PLYADD (K, A,B ,C) 

Ihis  subroutine  computes  C = A+B,  A,B  and  C are  Kth  order 
polynomials . 

SUBROUTINE  PLYSUB (K, A, B ,C) 

This  subroutine  computes  C = A-B,  A,  B,  and  C are  Kth  order 
polynomials . 

SUBROUTINE  IPYMUL (K1 , K2 , A, B , C) 

This  subroutine  computes  C = A x b A is  Kith  order  poly- 
nomial, B is  K2th  order  polynomial  and  C is  (Kl+K2-l)th  order 
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polynomial . 

SUBROUTINE  PLYDRV (K, A, B) 

Ihis  subroutine  computes  B = A' , A is  Kth  order  polynomial. 
2.4.4  Problem  Dependent  Subroutines. 

There  are  five  problem  dependent,  user  supplied  subroutines. 

Two  of  the  subroutines  compute  with  interval  arithmetic  and  three 
with  double  precision. 

SUBROUTINE  FUNC (XI , SOL, FSOL, L, K, KF) 

This  subroutine  H is  given  the  polynomial  representation  of 
XQ(t)  on  the  interval  [xl  (1)  ,xl  (2)  ] cind  it  confutes  the  poly- 
nomial representation  of  f(t,XQ(t)).  The  solution's  polynomials 
are  Kth  order  and  are  interval  valued. 

The  parameters: 

xl : - Interval  valued  vector  of  dimension  2.  It  gives  the  end 
points  of  the  interval  on  which  the  polynomial  representation  is 
valid. 

SOL:-  KxL  interval  array  with  the  coefficients  of  the  poly- 
nomials . 

FSOL:-  KFxL  array  with  the  coefficients  of  f(t,XQ(t)). 

L:-  The  size  of  the  system. 

KF:-  The  first  dimension  of  FSOL. 

SUBROUTINE  SETA (K ,L, XI , SOL, A, KA) 

This  subroutine  is  given  the  polynomial  representation  of  the 

solution  and  it  computes  the  polynomial  representation  of  A(t)  = 

f (t,x  (t)). 

X 0 
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The  uarameters: 


K,L,XI,SOL:-  The  same  as  in  FUNC 

A:-  A KAxLxL  interval  array  with  the  coefficients  of  A(t) . 
Remarks . 

1)  The  above  subroutines  should  be  translated  with  AUGMENT. 

2)  The  user  could  use  the  polynomial  manipulation  subroutines 
to  compute  f(t,XQ(t))  and  A(t)  . 

3)  Note  that  the  polynomials  are  of  the  form 

K 

P(T)  = y D(J)*(T-XI{1))**(J-1) 

J=1 

Therefore, the  polynomial  representation  of  the  independent  variable 
t on  the  interval  [XI  (1),  XI  (2)]  is  (XI  (1) ,1 ,0, . . . ,0) . 
SUBROUTINE  DEFUNC (T, Y, DY , STOR, IFLAG) 

This  subroutine  is  given  the  value  of  the  solution  at  T in 
Y and  it  confutes  f(T,Y)  and  stores  it  in  DY.  This  subroutine 
is  used  in  conputing  the  piecewise  polynomial  representation  of  the 
solution . 

The  parameters: 

T:-  Point  of  evaluation, 

Y:-  L dimensional  do\ible  precision  vector. 

The  value  of  the  solution  at  T . 

DY:-  L dimensional  double  precision  vector  contains  the 
values  of  f (T, Y) . 

STOR:-  Not  used. 

IFLAG:-  Not  xjsed. 
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SUBROUTINE  DERIVS (T , Y , DY , STOR, IFLAG) 

TTiis  subroutine  is  used  in  order  to  compute  the  fundamental 
solution.  It  is  given  in  Y the  values  of  one  column  of  Y(t) 
and  in  DY  the  values  of  A{t)Y. 

The  para^eters: 

T,Y,DY:-  The  values  of  Y(I),  I = 1,...,L  have  to  be  put  in 
STOR  (See  [49]  for  saving  options). 

IFLAG:-  not  used. 

SUBROUTINE  TRDERV (T,Y,DY, STOR, IFLAG) 

-1  T 

This  subroutine  is  used  in  order  to  compute  (Y  (t) ) . It  is 

-1  T 

given  in  Y the  values  of  one  column  of  (Y  (t))  and  stores  in 
T 

DY  the  vector  -A  (t)Y. 

The  parameters  are  the  same  as  in  DERIVS. 

-1  T 

The  differential  eqioations  that  Y(t)  and  (Y  (t) ) satisfy, 
depend  on  the  solution  x^Ct).  Therefore,  the  values  of  the  solu- 
tion at  the  knot  sequence  are  put  into  an  array  in  the  common 
bloclt  /SPLCFF/LXI ,C  ( • ) - LXI  is  the  number  of  points  the  solution 
is  computed  at  and  C is  LXIxL  double  precision  array  with  the 
solution  values  at  the  Icnot  sequence . 

In  order  to  facilitate  the  computation  of  point 

a subroutine  is  provided: 

SUBROUTINE  PPVLUl (C, LXI , YHT, X) 

The  parameters : 

LXI:-  Number  of  points  in  the  )cnot  sequence. 

C:-  LXIxL  array  with  the  values  of  the  solution  at  the  )<not 
sequence . 


-50- 


L:- 


The  size  of  the  system. 


YHT:-  L dimensional  array  with  the  values  of  the  solution 
at  X. 

X:-  Point  of  evaluation. 

2.4.5  The  Subroutine  APOSTR. 

The  main  subroutine  is  APOSTR.  It's  calling  sequence  is: 
SUBROUTINE  APOSTR (KA, KF , L , NPT, SOLTN, FUNC, SETA , B1 , B2 , NRMINV, ETANRM , 
WRKSPC) . 


The  parameters : 

KA:-  The  order  of  the  polynomials  in  the  matrix  A(t)  = 

f (t,x  (t)). 

X 0 

KF:-  The  order  of  the  polynomials  in  f(t,x^(t)) 

(Remember  that  x^ (t)  is  6th  order  piecewise  poly- 
nomial function) . 

L:-  The  size  of  the  system  of  differential  equations. 

NPT:-  The  number  of  points  in  the  knot  sequence. 

SOLTN:-  NPTxL  double  precision  array  with  the  values  of  the 
solution  knot  sequence. 

FUNC:-  Problem  dependent  subroutine. 

This  subroutine  is  given  the  polynomial  representation  of 
x^Ct)  and  it  computes  the  polynomials  f(t,x^(t))  for  calling 
sequence.  See  section  (2.4.4). 

The  subroutine  FUNC  should  be  declared  external  in  the  calling 
program. 
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SETA:-  Problem  dependent  subroutine. 

This  subroutine  is  given  the  polynomial  representation 
of  it  computes  the  matrix  polynomial  A(t)  = 

f^(t,Xp(t)).  For  calling  sequence  see  (2.4.4).  This  sub- 
routine should  be  declared  EXTERNAL  in  th  calling  program. 
B1,B2:-  LxL  double  precision  arrays  with  the  coefficients 

B,  and  B., . 

1 2 

NRMINV:-  Double  precision  value  of  the  bound  on 
ETANRM:-  Double  precision  value  of  the  bound  on 
WRKSPC:-  Double  precision  work  space  array. 

The  dimension  of  WRKSPC  should  be  at  least 
2*NPT*(L*L+2)+4*L*L*(38+3*KA) . 

Please  note 

The  knot  sequence  is  not  part  of  the  calling  sequence.  It 
must  be  put  into  the  common  block  /DPSTEP/  STEPS{-)  prior  to  the 
call  to  APOSTR. 
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2 . 5 Examples . 

In  this  section  we  give  two  examples  of  computational  existence 
proofs.  Analytical  existence  proofs  for  these  problems  are  not 
known. 


2.5.1  Example  1 . 

Consider  the  following  two-point  boundary  value  problem: 

ey"  = (y^  - (t-l)^)y'  0 1 t 1 1 (El) 

y(0)  = A , y(l)  = B . 


The  above  problem  was  suggested  and  analyzed  by  Howes  and 
Barter  as  a model  problem  for  nonlinear  problems  having  a contin- 
uous locus  of  singular  points  (see  [17]).  They  studied  the  asymp- 
totic behaviour  of  solutions  to  this  problem  as  e ->■  o"*^.  It  is 
not  hard  to  show  that  the  above  problem  has  at  least  one  solution 
for  any  value  of  e > 0.  However  the  existence  of  multiple  solu- 
tions was  not  ruled  out.  Moreover , Howes  and  Barter  showed  that  if 


0 < B r ^ ^ there  are  at  most  three  limit  solutions 


+ 1 
as  e->-0;y=A,  y=B  and  y 


/3 


Indeed  Francis  Sutton  [36] 
computed,  using  finite  differences,  three  solutions  for  e > 0 
small,  thus  implying  that  all  three  possible  limit  solutions  are  in 
fact  obtained. 

We  also  have  found  three  distinct  numerical  solutions  in  that 
range.  Our  numerical  solutions  were  obtained  by  a collocation 
method.  The  subroutine  LOBATO  by  deBoor  and  Weiss  [4]  was  used  to 
obtain  the  numerical  solutions.  Taking  A = .96,  B = .001  amd 
E = 1/15  the  above  aposteriori  error  analysis  was  used  to  estab- 
lish the  existence  of  (at  least)  three  distinct  solutions. 
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Moreover,  guaranteed  error  bounds  were  obtained. 

Of  course,  the  theory  does  not  rule  out  the  existence  of  more 
solutions  for  e >0.  However,  the  theoretical  results  of  Howes 
and  Parter  combined  with  Sutton's  computational  results  and  our 
existence  re;;ults  for  e = 1/15  seem  to  support  Sutton's  conjecture 
that  Equation  (El)  has  exactly  3 solutions  for  all  0 < e < , for 

some  Cq  > 0 . 

Although  e = 1/15  is  not  very  small  Equation  El  is  already 
"stiff".  As  a matter  of  fact  solution  3 (see  figure  2)  was 
"stiffer"  than  solutions  1 and  2.  In  order  to  get  reasonable  error 
bounds  we  had  to  modify  the  original  programs  so  that  it  would  be 
possible  to  subdivide  the  interval  into  unequal  subintervals . The 
bounds  on  the  residual  were  reduced  by  4-5  orders  of  magnitude  by 
putting  more  points  in  the  interval  where  the  function  and  it's 
derivative  change  very  fast.  Solutions  1 and  2 were  computed  using 
201  points  and  solution  3 lising  301  points.  The  bounds  we  obtained 
are  given  in  Teible  1.  We  now  turn  to  our  computations: 

Let  us  rewrite  equation  (El)  as  a first  order  system  with 
R = 1/e. 


the  matrix  of  partial  derivatives  is  given  by: 
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A(t)  = 


R(yJ  - (t-1)^) 


The  second  partial  derivatives  are: 


3yi9y2 


and 

there  fore 


!!!2 

•^2 

3 y-, 


a'f. 


= 0 


= 2Ry, 


3"y2  = 


= 2ry, 


sup 


F"  (x) 


^ 2R(2j’  |y  (s)  |ds  + / |y  (s)  |ds)  + 6 R r . 

0 0 ^ 


Therefore,  the  bound  on  the  second  Frechet  derivative  can  be  com- 
puted by  the  above  formula.  Since  y^  and  y^  do  not  change  sign 
on  [0,1]  the  above  integrals  can  be  coirputed  exactly  (Remember 
that  y^  and  y^  are  piecewise  polynomial  functions . ) . 

In  figure  2 we  have  plotted  the  three  different  solutions. 

The  error  bounds  we  obtained  are : 


K 

n 

B 

h 

^0 

77.5179 

2 .40683x10”^ 

4.4769x10^ 

8.3526 xio"^ 

2 .4069xlo”^ 

60.9082 

1.24 32 3x10"^ 

5.0692x10^ 

3.8386 xlo“^ 

1.24326xlo"® 

37.98871 

5. 39276xlo“^ 

1.7473x10^ 

3.5795x10"'^ 

5. 3937x10“^ 

Table  1 
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Figure  2 


2.5.2  Example  2 . 

This  exaitple  was  mentioned  in  the  introduction.  These  differ- 
ential equations  describe  the  flow  between  two  parallel  infinite 
disks  rotating  about  a common  axis.  The  equations  are: 
h^'^  + h h'"  + g g'  =0 


(E2) 


g"+hg'  -h'g=0 
h (0)  = h'  (0)  = h(l)  = h'  (1)  = 0 
g(0)  = , g(l)  = . 


Althou^  the  above  problem  has  attracted  considerable  attention 
(See  McLeod  [22]  and  the  references  there),  existence  proofs  for 
solutions  outside  a small  set  of  values  of  are  not  known 

(see  Elcrat  [14]  and  McLeod  [22]). 
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Hastings  proved  existence  provided  and  are  suffi- 

ciently small. 

Elcrat,  using  a fixed  point  argument  proved  that:  If 

a)  e and  0 ^ 

or 

2 2 2 

b)  e (0,?^^]  and  ^g  “ ^2  ' 

45  4 -1/2 

where  C = — (2)‘  (exp(.25)  + exp(— ) - 1)  = 1.5  then  equation 

(E2)  has  a solution. 

McLeod  cind  Parter  [28]  proved  existence  of  solutions  for  the 
counter  rotating  case  They  proved  existence  of  solu- 
tions for  all  values  of  > 0,  Proof  of  existence  in  all  other 

cases  is  not  known. 

Many  people  have  conputed  solutions  outside  the  range  where 
the  existence  of  solutions  is  known.  Moreover,  multiple  numerical 
solutions  were  obtained.  However,  none  of  the  computational  papers 
give  any  error  analysis.  Although  many  of  these  computations  prob- 
cibly  give  reasonable  approximation  to  solutions  of  equations  (E2)  , 
there  is  at  least  one  case  where  a computed  solution  was  proven  to 
be  incorrect.  (See  McLeod  and  Parter  [28].) 

By  using  an  a posteriori  error  analysis  , the  existence  of  a solu- 
tion, outside  the  range  where  the  previously  known  results  apply, 
has  been  guaranteed.  Although  we  have  not  proven  the  existence  of 
multiple  solutions  (when  and  are  large  enough),  we  believe 

that  with  enough  time  and  computing  power  (money)  one  could  use  our 
method  to  establish  the  existence  of  multiple  solutions.  As  a by- 
product of  Elcrat 's  proof  the  solutions  he  obtains  are  "monotone" 
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in  the  sense  that  g'  is  of  one  sign.  On  the  other  hanci,  McLeod 
and  Parter  showed  t^.at  if  there  is  a solution  in  the  case  where 
and  are:  positive,  large  enough  and  far  apart,  g*  of  that 

solution  must  change  sign.  Moreover,  numerical  solutions  where  g' 
does  change  sign  were  obtained  (for  example  see  Cerutti  [6]).  It 
will  be  interesting  to  prove  the  existence  of  such  solutions. 

Note  thar  the  above  problem  is  "stiff"  and  it  gets  "stiffer" 
as  I I and  | | grow . 

This  numerical  solution  was  computed  with  subroutine  LOBATO  at 

101  equally  spaced  points.  We  computed  with  (2^  = 7 and  (2^^  = 1 . 

2 2 

Note  that  (2„  - (2,  = 48. 

0 1 

The  bounds  we  obtained  are: 


K 

n 

B 

h 

^0 

4 

2.0154x10“^ 

1. 60437x10 

.1293406 

2.16597xlo“® 

In  Figure  3 we  have  plotted  the  values  of  and  (2^^  for  which 

existence  proofs  are  known. 
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2 .6  Conclusions,  Remarks,  etc. 


In  this  work  we  have  demonstrated  a way  to  compute  aposteriori 
error  bounds  for  polynomials  two-point  boundary  value  problems. 

Some  of  the  difficulties  of  other  methods  were  overcome.  However, 
the  problem  is  by  no  means  completely  solved.  We  list  below  prob- 
lems and  questions  left  to  be  answered: 

a)  The  bound  on  h = |i  II  use,  is  not  the  best  possible 

one.  We  suggested  two  other  methods  for  computing  n . These 
methods  have  to  be  further  investigated.  Both  methods  suffer  from 
some  difficulties.  The  first  method  uses  the  explicit  form  of  the 
Green's  function  and  requires  the  manipulation  of  high  degree  poly- 
nomials. The  second  method  looks  more  promising,  however,  a large 
amount  of  work  may  be  necessary  in  order  to  make  the  residual,  that 

is  x^(t)  - /^[A(s)x^(s)  - r(s))ds  very  small.  In  any  case  more 
a 

experiments  are  needed  before  one  could  obtain  better  methods  for 
aposteriori  error  bounds. 

b)  We  have  chosen  to  interpolate  the  solution  x^(t),  y(t)  and 

Y ^(t)  by  6th  order  polynomials.  However,  there  are  no  compelling 
reasons  for  doing  so.  More  experiments  are  needed  to  obtain  some 
ideas  and  to  gain  some  insight  into  good  strategies  for  choosing 
..  erpolation  schemes. 

c)  The  implementation  of  interval  arithmetic,  we  use,  is  very  slow. 
A necessary  condition  for  making  interval  arithmetic  practical  for 
routine  use  is  having  the  arithemtic  done  by  hardware  and  not  by 
software.  The  prospects,  in  the  near  future,  of  having  interval 
arithmetic  as  a standard  hardware  option  are  not  good.  However, 


-59- 


a new  generation  of  microprogramable  computers  may  make  the  use  of 
interval  arithmetic  practical. 

d)  We  use  interval  arithmetic  in  order  to  take  round-off  error  into 
account.  If  one  wishes  to  prove  the  existence  of  solutions  one  is 
forced  to  take  round-off  error  into  account.  However,  it  is  our 
experience  that  computations  with  double  precision  arithmetic  give 
the  same  results  and  round-off  errors  have  negligible  effects.  Thus, 
our  method  could  be  used  to  compute  reliable  eri-'r  bounds  although 
the  results  would  not  be  completely  rigorous  existence  proof. 

e)  The  size  of  the  residual  or  the  size  of  ||  II  ori  each 

subinterval  could  be  used  to  decide  how  to  redistribute  or  refine 
the  subdivision  points.  More  experiments  are  needed  in  order  to 
find  strategies  based  on  the  above  information.  Also  one  should 
compare  these  methods  with  other  methods  suggested  in  the  literature 
(see  [2],  [7],  [45],  [48]). 

f)  As  was  said  in  the  introduction  our  method  could  be  extended  to 

problems  that  are  not  polynomials.  If  the  function  f(t,y)  of 
equation  (1)  is  a factorable  function  of  t and  y (see  Part  A of 
our  thesis),  then,  since  ® piecewise  polynomial  function, 

f(t,x„(t))  and  A(t)  = f (t,x„(t))  are  piecewise  factorable  func- 

0 X 0 

tions . Since  one  knows  the  Taylor  series  expansion  of  Xg(t)  on 
each  subinterval,  one  can  compute  Taylor  series  approximations  to 
f(t,XQ(t))  and  A(t)  . 

If  f(t)  is  the  approximation  to  f(t,x^(t))  and  A(t)  is 
the  approximation  to  f^(t,XQ{t)),  then: 


I 
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i)  Using  interval  techniques  (see  Moore  [30]  ch.  10)  one  can 

bound  ||  f (t) -f  (t  ,x^  (t)  )||  and  ||A(t)  - A(t)||  . Also  if  one  ta)tes 

enough  terms  in  the  series  expansion  and  the  subintervals  are  small 

enough,  these  bounds  will  be  small. 

ii)  Since  f(t)  is  a piecewise  polynomial  function  we  can  use 

it  in  order  to  compute  an  approximation  to  r(t)  . 

Recall  that  r(t)  = X|^(t)  - ^^(a)  - f (s  ,x^  (s) ) ds . If  r(t)  = 
t ~ a _ 

x^Ct)  - XpCa)  - / f(s)ds  then  r(t)  - r(t)  = -/  if  (s , x^(.s) ) - 
~ a a 

f(s))ds  therefore 

II  1^11  1 II  J^ll  + (b-a)  • II  f (tjX^Ct)  )-f  (t)  II  . 

iii)  Similarly,  A(t)  can  be  used  to  compute  approximation  to 

X,  (t)  and  (t)  . 

1 2 

Recall  that 


X^(t)  = V(t)  - V(a)  - A(s)V(s)ds  . 

a 

If 

Xj^(t)  = V(t)  - V(a)  - A(s)V(s)ds 

a 

then 

II  ^-xjl  1 I|a-a||  • (/^  |v{s)  Ids)  . 

a 

The  same  way,  if  ^2^^^  “ W(t)  - W(a)  + / W(s)A(s)ds  then 

II  • 

a 

In  turn  the  bounds  on  ||  X^^  ||  cind  ||  X^  ||  can  be  used  to  compute 
bounds  on  ||  E^||  and  ||  ||  . (See  section  (1.3.2)). 

Therefore,  an  algorithm,  similar  to  the  one  used  for  polynomial 
equations,  can  be  used  for  aposteriori  error  analysis  of  factorable 
two-point  boundary  value  problems. 

g)  It  is  not  hard  to  see  how  one  can  extend  our  method  for  problems 
with  nonlinear  boundary  conditions. 
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